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ABSTRACT 


The  report  contains  the  numerical  Implementation  of  a theo- 
retical algorithm  that  generalizes  the  Euler  equations  for  Irregular 
flows.  The  main  features  of  the  algorithm  are  described,  and  the 
relevant  equations  are  summarized.  The  spatial  region  occupied  by 
the  fluid  Is  then  discretized,  and  the  numerical  quadrature  of  the 
theoretical  algorithm's  governing  equations  Is  effected.  Computa- 
tional results  are  given  for  the  following  hydrodynamic  free-boundary 
problems:  (a)  the  fall  from  rest  of  a liquid  with  a profoundly  non- 
linear Initial  free  surface;  (b)  the  motion  from  rest  of  a liquid 
whose  Initial  free  surface  Is  only  slightly  removed  from  Its  equi- 
librium position;  and  (c)  the  collisions  of  streams  of  fluid  with 
formation  of  a jet.  A program  listing  Is  given  In  the  appendix. 
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1.  INTRODUCTION 


As  part  of  a continuing  investigation,  we  have  reformulated 
inviscid  hydrodynamics  in  a manner  that  is  natural  for  the  study  of 
the  free-surface  probl'*".  The  ultimate  purpose  of  the  investigation 
is  to  compute  in  an  effii  >t4.at  and  reliable  way  the  phenomena  at- 
tendant to  the  motion  of  a rigid  body  in  water  with  a free  surface. 

The  considerations  that  have  guided  our  reformulation  of 
hydrodynamics  and  the  reduction  of  our  generalized  theory  to  the 
classical  theory  when  the  flow  variables  are  sufficiently  smooth 
were  discussed  in  Ref.  1,  and  they  need  not  be  repeated  here.  The 
purpose  of  this  report  is  to  see  directly  what  the  implications  of 
the  generalized  hydrodynamics  are  for  th^  numerical  solution  of 
problems  with  hydrodynamic  free  surfaces  and  to  present  our  numer- 
ical results. 
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Nevertheless,  some  recapitulation  of  the  essential  ideas  of 
the  theory  is  in  order.  Hydrodynamics  as  reformulated  does  not 
take  on  quite  a Lagranglan  or  an  Eulerlan  mode  but  has  some  of  the 
advantages  of  both.  The  advantage  our  formulation  shares  with  the 
Eulerlan  one  is  that  the  equations  have  as  independent  variables 
the  ones  of  direct  physical  significance  — space  and  time.  The 
advantage  it  has  in  common  with  the  Lagranglan  approach  is  that 
time-dependent  free-surface  problems  may  be  solved  by  solving  a 
system  of  equations  on  a domain  that  is  Independent  of  the  time. 

The  evolution  of  an  incompressible  flow  is  seen  as  the 
solution  of  a set  of  hyperbolic  conservation  laws  subject  to  a 
constraint.  The  hyperbolic  conservation  laws  are  Just  the  equa- 
tions of  mass  and  momentum  conservation.  In  the  absence  of  the 
constraint,  they  are  the  equations  of  a perfectly  compressible 
(pressureless)  fluid.  (In  N dimensions,  we  also  refer  to  these 
as  the  N-dlmenslonal  inviscid  Burgers  equation.)  The  constraint 
is  a one-sided  constraint  on  the  density  that  expresses  incompress- 
ibility by  establishing  an  upper  bound  on  the  amount  of  fluid  that 
may  occupy  any  volume  of  space. 

In  practice,  in  our  theory  the  solution  of  the  combined 
conservation  laws /constraint  is  achieved  by  going  from  one  time 
to  a slightly  enhanced  one  in  a "split-step"  scheme,  in  which  at 

Ref.  1. J.  C.  W.  Rogers,  "Incompressible  Flows  as  a System 
of  Conse|^atlon  La%ra  with  a Constraint,"  Semlnaires  IRIA,  Analyse 
et  Controle  de  Systraes,  1978. 
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first  the  conservation  laws  are  solved  as  If  there  were  no  constraint, 
and  then  the  constraint  Is  satisfied  In  a manner  that  retains  the  con- 
servation of  mass  and  momentum.  Thus,  If  the  solution  of  the  evolu- 
tionary problem  at  a given  time  Is  thought  of  as  obtained  through  the 
action  of  a nonlinear  semigroup  on  the  Initial  data,  our  theory  pro- 
vides an  approximation  to  Che  semigroup. 

In  Section  2,  we  will  summarize  the  governing  equations  that 
are  solved  at  each  time.  Some  numerical  quadratures  that  enable 
these  equations  to  be  put  In  a flnltary  form  suitable  for  computer 
manipulation  are  given  In  Section  3.  Section  4 gives  numerical  re- 
sults for  three  sample  problems  that  have  been  run.  The  first  of 
these  follows  the  time  evolution  of  an  Incompressible  fluid  that 
was  Initially  at  rest  and  whose  Initial  free  surface  was  distorted 
In  a profoundly  nonlinear  manner.  In  the  second  problem,  a com- 
parison Is  made  between  the  computer  results  and  the  predictions 
of  linearized  water  wave  theory.  The  third  example  Illustrates 
the  computation  of  phenomena  associated  with  the  collision  of  two 
streams  and  the  formation  of  a jet.  Here  there  Is  no  exact  theory 
to  compare.  Instead,  the  results  of  numerical  computations  are 
compared  for  different  mesh  sizes  and  time  steps.  All  the  numerical 
examples  presented  here  are  for  problems  with  two  Independent  space 
variables.  Section  5 discusses  some  of  the  limitations  of  our 
method,  possible  Improvements,  and  some  of  our  plans  for  further 
computational  work.  In  the  belief  that  nothing  removes  ambiguity 
like  an  explicit  statement  of  the  steps  we  have  gone  through  to 
Implement  our  theoretical  Ideas  numerically,  we  have  Included  a 
program  listing  In  Appendix  A. 

Several  observations  are  In  order.  First,  the  generalized 
hydrodynamic  theory  Is  by  no  means  completed.  Since  our  theory  Is 
given  In  constructive  fashion  through  an  algorithm  to  determine  the 
flow  at  any  time  In  terms  of  the  flow  a time  step  earlier,  the  theory 
will  be  acceptable  only  when  we  have  obtained  definitive  results  for 
the  construction  as  the  time  step  goes  to  zero.  Such  results  will 
have  to  Include  the  regularity  of  the  flow,  convergence  of  the  con- 
struction to  a semigroup  In  the  appropriate  function  spaces,  and 
existence  of  the  flow  globally  In  time.  Indications  to  date  are 
that  the  conclusion  of  these  tasks  In  a satisfactory  manner  will  be 
concomitant  with  the  development  of  a theory  of  Invlscld  hydrodynamic 
turbulence  (Ref.  2). 

In  spite  of  this  reservation,  the  conasencement  of  calculations 
with  the  generalized  theory  Is  by  no  means  premature.  For  one  thing, 
the  fact  that  the  general  formulation  reduces  to  the  usual  one  when 

Ref . 2 . J.  C.  W.  Rogers,  "Stability,  Energy  Conservation, 
and  Turbulence  for  Water  Waves,"  Semlnalres  IRIA,  Analyse  et  Controle 
de  Systemes,  1978. 
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the  flow  variables  are  sufficiently  saooth  guarantees  that  the  status 
of  the  fundamental  questions  i;loned  above  Is  no  worse  In  our 
theory  than  In  classical  hydrodynamics.  Accordingly,  calculations 
based  on  our  approach  are  no  less  sound  than  those  based  on  conven- 
tional approaches.  To  the  contrary,  the  fact  that  the  general  theory 
makes  sense  In  e wider  verlety  of  circumstances  and  under  less  regu- 
larity requirements  on  the  flow  than  the  usual  one  suggests  that  per- 
haps a greater  presumption  of  success  In  resolving  the  basic  questions 
of  existence,  regularity,  and  turbulence  may  be  attached  to  the  gen- 
eral theory  than  to  the  classical  one. 


A second  observation  relates  to  the  fact  that  In  most  cases 
of  interest  the  flow  of  an  unstratlfled.  Incompressible,  Invlscld 
fluid  Is  Irrotatlonal.  The  general  theory  applies  to  rotational  as 
well  as  Irrotatlonal  flows.  However,  for  the  Important  special  case 
of  Irrotatlonal  flows,  significant  savings  In  computational  time  may 
be  effected  In  the  classical  theory  through  the  use  of  Integral  equa- 
tions. An  unfinished  task  Is  to  examine  the  ramifications  of  the 
general  theory  for  Initially  Irrotatlonal  flows  and  to  develop,  to 
whatever  extent  It  Is  possible,  variations  on  our  algorithm  that  re- 
tain its  essential  character  but  effect  the  computational  savings 
anticipated  for  this  special  case. 


Therefore,  we  do  not  consider  the  algorithm  we  have  studied 
to  be  final,  by  any  means,  for  the  majority  of  cases  of  practical 
Interest.  Accordingly,  we  have  not  expended  great  energy  In  trying 
to  perfect  numerically  the  algorithm  In  the  form  in  which  it  now 
stands.  From  a computational  point  of  view,  this  report  should  be 
viewed  as  a preliminary  report  whose  purpose  is  to  show  the  essen- 
tial correctness  of  our  approach.  In  practice  as  well  as  In  theory. 

We  hope  these  coomants  put  the  accompanying  numerical  work  In 
a proper  perspective.  While  we  have  not  been  deliberately  careless 
In  carrying  out  the  numerical  quadratures,  we  have  also  not  made  a 
number  of  refinements  that  could  have  been  made.  (These  points  are 
discussed  more  fully  In  Section  5.)  Thus,  to  some  extent,  the  theory 
has  been  tested  under  rather  adverse  conditions.  Our  expectation  Is 
that.  If  the  numerical  Implementation  of  the  theory  thus  handicapped 
yields  at  all  reasonable  results,  more  can  be  expected  of  a similar 
approach,  executed  with  greater  care.  We  have  ventured  forth  In 
this  manner,  guided  by  our  faith  that  the  theory,  free  as  It  Is  of 
the  comparatively  severe  regularity  requirements  of  the  classical 
hydrodynamic  theory,  has  an  essential  robustness  that  enables  it  to 
carry  the  weight  of  even  a crude  numerical  quadrature. 
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2.  GOVERNING  EQUATIONS 


Let  D be  a domain  occupied  by  the  fluid.  In  the  problems 
treated  In  this  paper,  D Is  Independent  of  time,  and  the  boundary 
dD  Is  rigid.  D need  not  be  bounded,  p will  denote  the  fluid  den- 
sity. The  density  constraint  is 

P S Pq,  (1) 

where  p^  Is  the  density  of  the  fluid  In  Its  Incompressible  (liquid) 
phase. 

The  velocity  field  Is  denoted  by  u(x,t),  and,  of  course,  the 
momentum  density  Is  p(x,t)u(x,t) . The  evolutionary  problem  Is  to 
find  p(x,t)  and  p(x,t)u(x,t)  given 

p*^(x)  - p(x,0)  (2a) 

and 


p®(x)u®(x)  - p(x,0)u(x,0). 


(2b) 


We  write  the  solution  symbolically  as 

[p(x,t),  p(x,t)u(x,t)]  - S(t)  (p®,  p^u®),  (3) 

where  S(t)  Is  a nonlinear  semigroup  satisfying 

S(ti  + tp  - S(tj^)  S(.t^).  (4) 

(In  some  turbulent  flows,  the  system  may  evolve  stochastically  even 
though  the  Initial  state  Is  uniquely  prescribed,  and  In  that  case 
Eq.  4 Is  a statement  of  the  Markov  property  of  the  time  evolution 
of  the  flow.)  A particular  case  of  Eq.  4 Is 

S(t)  - [S(t)]",  (5) 

where  t ■ nt  and  x Is  called  the  time  step. 
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The  algorithm  we  have  Introduced  (Ref.  1)  gives  an  approxl' 

matlon  to  S(t),  which  we  denote  by  S(t).  S(t)  Is  determined  as 
follows:  we  suppose  we  are  given  p(x)  and  p(x)u(x),  with  p(x) 
satisfying  the  constraint  L4.  1.  First,  we  solve  the  hyperbolic 
"conservation"  laws. 


with  Initial  conditions 


C(x,0)  ■ p(x),  xcD,  C(x,0)n(x,0)  - p(x)u(x),  xeD 


and  boundary  conditions 


where  n Is  the  unit  outward  normal  to  dD.  In  terms  of  the  solution 
of  Eqs.  6,  7,  and  8 at  t t,  we  define  the  qtiantltles 


This  completes  the  first  part  of  our  "split-step"  scheme 


Next,  we  solve  the  one-phase  Stefan  problem 


subject  to  the  Initial  condition 


r f 

t 

it 
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and  the  boundary  condition 

V8*n  - 0,  (x,a)c  3D  x (0,»). 


(lOd) 


As  a-x»,  6(x,a)  approaches  a steady-state  value  that  we  denote  by 
p(x);  that  Is, 


P(x)  H 0(x,a),  xeD. 


Note  that  p(x)  satisfies  Eq.  1. 


Continuing,  we  define 


v(x)  = y*  f [0(x,a)]  do,  xeD, 


(11) 


(12) 


where  0 Is  given  by  Eq.  10,  and  we  determine  u(x)  as  the  solution 
of 


I 


and  subject  to  the  boundary  conditions 


With  p(x)  and  p(x)u(x)  thus  obtained,  we  define  S(t)  as  the  operator 
auch  that 


The  solution  of  the  Stefan  problem  (Eq.  10)  and  the  linear 
elliptic  equation  (Eq.  13)  Is  straightforward.  It  is  the  solution 
of  the  Stefan  problem  that  determines  the  time  development  of  the 
free  boundary.  Explicit  numerical  algorlthaw  to  solve  these  prob- 
lems follow. 
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For  the  Stefan  problem,  we  use  a variation  on  an  algorithm 
for  the  general  problem 

+ Lf(u)  - 0,  u(x,0)  ••  Uq(x),  (15) 

when  the  semigroup  associated  with  L, 

S(t)  - e^^,  (16) 

1 00 

Is  contractive  In  L and  L (Ref.  3).  The  algorithm  approximates 
u(x,nh)  by  u°(x),  obtained  by  making  the  substitutions 


„ . u-C.)  - 

t h 


(17a) 


and 


, . S(Ah)  - 1 

^ - Ah '*• 


(17b) 


where  A la  a positive  constant,  In  Eq.  15.  Thus, 


u”'*'^(x)  - u*‘(x)  - f[u“(x)]  S(hA)  f[u“(x)] 

and 

u°(x)  - Uq(x).  (18) 

1 00 

The  algorithm  (Eq.  18)  Is  stable  In  L and  L If  f satisfies 

0 « f(u)  - f(v)  » A(u  - v)  for  u - V ^ 0.  (19) 

Kaklng  the  obvious  substitutions  and  noting  the  boundary  condition 
(Eq.  lOd),  we  obtain  an  algorithm  to  solve  the  Stefan  problem  (Eq. 
10)  by  setting  A - 1: 
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0“(x)  P'  0(x,nAa), 


(20a) 


0°(x)  - p(x). 


(20b) 


0“'^^(x)  - 0“(x)  - f[0“(x)]  + S(Aa)  fl0“(x)], 


(20c) 


where  ^(x)  is  obtained  from  0*'(x)  by  reflecting  values  of  0*'(x) 
symmetrically  across  the  boundary  9D: 


0®(x)  - 0 (x),  xeD, 


(20d) 


0”(x+ne)  ■ 0*'(x-ne),  e ^ 0,  xc9D. 


(20e) 


In  this  case,  L « -A  and 


tS(h)u](x)  - 


1 r -(x-x')' 

zjn  I * 


u(x’)dx’. 


(4Trh)"'‘  -^N 


There  are  many  ways  to  solve  the  linear  elliptic  boundary 

value  problem  (Eq.  13)  for  p u.  In  terms  of  the  solution  of  the 
parabolic  problem 


* - A(f-  i|<), 


(22a) 


♦(x.O)  - p(x)u(x) 


(22b) 


we  get 


.u(.)  ./ 


♦(x,Y)e"'dY. 
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An  approximate  solution  of  Eq.  22  Is  obtained  by  making  substitu- 
tions In  Eq.  22a  of  the  sort  Indicated  In  Eq.  17.  Thus,  If 


(x)  fe  <l)  (x.nAy) , 


(24a) 


we  get 


<|»°(x)  - p(x>u(x) 


iw. 


(24b) 


and 


'l»”^^(x)  - ♦“(x)  - 7 ^ + 7 S(AAy)  (^  ijl®)  . (24c; 

A Pq  a Pq 

Here,  In  accordance  with  the  boundary  conditions  (Eq.  13b),  we  ob- 
tain V*  from  by  reflecting  the  components  of  parallel  to  the 

boundary  9D  symmetrically  and  the  component  of  <1^°  perpendicular  to 
3D  antisymmetrically: 


^"(x)  - ip^(x),  xeD, 


and 


(^  X n)(x+ne)  ■ (p  x n)(x-n€),  e a 0,  xe3D, 


(v  •n)(x+nc)  ■ -(v  •n)(x-ne),  e * 0,  xe3D. 


The  scheme  (Eq.  24c)  Is  stable  In  L and  L If 


v(x). 


Pq  xeD 


(24d; 


(24e: 

(24f; 


(25! 


S(AAy)  Is  given  by  Eq.  21. 
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The  function  v(x),  which  appears  In  Eq.  24c  and  Is  defined 
In  Eq.  12,  is  approximated  by 

00 

v(x)  « Ao  fie^Cx)].  (26) 

n-0 


2 

The  term  - In  Eq.  13a  had  Its  origin  In  the  momentum  associated 

with  the  redistribution  of  mass  upon  satisfying  the  constraint  (Eq. 
1)  (Ref.  1).  Thus,  this  term,  which  appears  In  Eq.  24b,  Is  computed 
approximately  from  a formula  that  reflects  Its  origin: 


- -^v  (x) 

T 


^ (4ffAa) 


n“0 


(27) 


In  contrast  to  the  solution  of  Eqs.  10  and  13,  the  solution 
of  the  hyperbolic  conservation  laws  (Eqs.  6 through  9)  poses  larger 
theoretical  problems,  as  it  Is  closely  connected  with  the  origins 
of  turbulence  and  In  the  general  case  requires  the  enlargement  of 
the  class  of  acceptable  solutions  to  stochastic  flows  In  order  to 
be  well  posed  (Ref.  2).  At  this  point  we  come  close  to  the  current 
limitations  of  our  hydrodynamic  theory.  In  particular,  we  do  not 
now  have  a reliable  algorithm  to  determine  the  evolution  of  the 
flow  in  probability.  Accordingly,  we  shall  assume  that  all  flows 
studied  In  this  report,  whether  turbulent  or  not,  evolve  determin- 
istically. 

The  purpose  In  regarding  Eqs.  6 through  9 as  a "conservation" 
law  Is  to  provide  a guide  for  determining  the  "weak"  solution  of 
Eq.  6 when  the  Initial  data  p(x)  and  p(x)u(x)  lack  sufficient  reg- 
ularity for  a classical  solution  to  exist  for  all  te(0,T).  An  ap- 
proximate solution  of  this  conservation  law  can  be  given  in  terms 
of  a distribution  function  F(x,v,t)  satisfying  the  equation 

+ vTF  - g 1^  - 0,  (x,v,t)eD  x R**  x(0,t),  (28) 

z 


the  Initial  condition 


F(x,v,0)  • p(x)  d(v-u(x)],  (x,v)eD  x R*’, 


(29) 
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and  the  boundary  condition 


F(x,v,t)  ■ F(x,v-2nvn,t) , (x,v,t)e3D  x R x (0,t). 
p(x)  and  p(x)u(x)  are  given  approximately  by 


(30) 


p(x)  - y*F(x,v,T)dv,  xeD,  p(x)ij(x)  - J"  F(x,v,t)v  dv,  xeD.  (31) 

Note  that  Eq.  28  Is  a linear  equation  whose  solution  can  be  written 
explicitly.  Because  of  the  boundary  condition  (Eq.  30),  the  char- 
acteristics of  the  equation  satisfy 


dt 


(32a) 


^ 2nvn  6(t-tj^), 

where  the  times  {t^}  are  the  times  when 
x(tj^)e3D 


(32b) 


(32c) 


and  n Is  the  outward  normal  to  3D  at  x(tj^).  The  use  of  the  distri- 
bution function  to  solve  the  conservation  law  (Eqs.  6 through  9) 
has  some  similarity  to  the  construction  of  solutions  of  another 
hyperbolic  conservation  law  through  the  superposition  of  solutions 
of  linear  equations  (Ref.  A). 


Ref.  A.  J.  C.  W.  Rogers,  "An  Algorithm  for  e Hyperbolic 
Free  Boundary  Problem,"  APL/JHU  TG  1309,  Kay  1977. 
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3.  NUMERICAL  SOLUTION  OF  THE  EQUATIONS 


Let  X and  z be  the  two  Independent  spatial  variables.  We 
shall  Indicate  how  the  equations  given  In  the  last  section  are 
treated  numerically  when 

D - {(x,z)  I z > 0,  0 < X < X}.  (33) 

We  approximate  D by  the  computational  domain  D : 

c 


- {(x,z)  I 0 < X < X,  0 < z < Z}, 


and  we  partition  Into  rectangles  by  a grid  of  lines 


TT<«- 


*1+1/2^  ^*”*j+l/2^  ■°*  OCIKI,  OKjAJ 


with 


*1/2  - 0*  *1-1/2  " *1+1/2  1 K 1 K I,  Xj^^/2  “ ** 

*1/2  * *j-l/2  *j+l/2  1 K j K J, 


(34) 


(35a) 


(35b) 


The  fundamental  dependent  variables  are  m 
rectangle 


Ij 


> 


the  mass  In  the 


®lj  = ^*1-1/2 » *1+1/2^  ^^J-1/2’  *J+l/2^’ 

And  (resp.  x-component  (reap,  z-component)  of  momen- 

tum In  the  same  rectangle.  Each  of  these  quantities  should  be 
labeled  by  another  Index  to  distinguish  the  time  at  which  It  Is 
measured.  However,  for  purposes  of  simplicity  and  In  the  spirit  of 
the  analytical  description  of  the  algorithm  In  Eqs.  6 through  14, 
we  do  not  carry  this  Index  along  and  only  Indicate  how  to  find  these 
quantities  at  one  time  from  their  values  one  time  step  beforehand. 
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We  first  show  the  effect  of  solving  the  hyperbolic  conser- 
vstlon  lew  (Eqs.  6 through  9)  on  the  quantities  ^xlj'  ^zlj * 

Here  we  use  the  approximate  form  of  the  algorithm  In  Eqs.  28 
through  31.  The  first  step  Is  to  take  account  of  the  effect  of 
gravity  on  u,..; 


^Ij  " “zlj  ■ 


Next  we  compute  the  velocities  u^^  and  w^^ . Given  an  appro- 


priate small  number  e>0,  we  have 


.y<e 


(38.) 


"ij" 


"IJ*' 


(38b) 


The  mass  density  Is 


O • " , 

‘’ij  Ax^AZj  * 

where 


Ax.  - X, 


1 *1+1/2  ■ *1-1/2’  " *j+l/2  " *j-l/2* 


Our  numerical  solution  of  Eqs.  28  through  31  proceeds  as  If 
p(x),  u(x),  and  w(x)  are  each  constant  In  rectangle  with  values 

^ij’  “ir  and  Wj^j,  respectively.  The  boundary  conditions  (Eq.  30) 
hold  on  the  rigid  part  of  9D: 


- 21  - 


j 


THE  X)HNS  HOPKINS  UNIVERSITY 

APPUEO  PHYSICS  LABORATORY 
LAUREL  Maryland 


r = {0}  X (0,z)  U (0,X)  X {0}  U {X}  X (0,Z).  (41) 

On  3D^-r,  we  allow  fluid  to  leave  the  computational  region  and  never 
return. 


With  these  boundary  conditions,  the  problem  becomes  determi- 
nate. In  Eq.  37,  we  have  taken  account  of  the  effect  of  gravity. 
There  remains  to  solve  an  equation  like  Eq.  28  with  g > 0 and  the 
boundary  conditions  (Eq.  30).  The  equations  of  the  characteristics 
are  Eqs.  32  with  g ■ 0.  For  our  computational  domain  the  solution 

of  these  equations  Is  straightforward.  One  need  only  translate  each 
point  of  each  rectangle  along  Its  appropriate  characteristic  for 

a time  t,  find  Its  new  location  In  the  computational  grid,  and.  In 
addition,  ascertain  the  number  of  reversals  of  the  normal  component 
of  velocity  that  have  taken  place  at  F,  In  accordance  with  Eqs.  32. 


Different  cases  arise,  according  to  whether  the  characteristic 
for  a point  of  a rectangle  has  or  has  not  been  reflected  at  z ■ 0, 

or  whether  the  characteristic  has  left  the  computational  grid.  A 
similar  situation  arises  with  regard  to  reflection  of  characteristics 
at  X ■ 0 and  x ■ X,  In  particular,  whether  the  total  number  of  such 
reflections  Is  even  or  odd. 


Let  xeD^  be  a point  In  the  original  computational  domain  and 
let  x*(x)eD^  be  Its  new  location  after  time  t.  If  Its  characteristic 

does  not  leave  D^.  We  may  map  each  such  point  x Into  a point  x(x) 

of  the  rectangle  (0,2X)  x (0,2Z)  according  to  the  following  prescrip- 
tion. 


X*(x) 

^1 

even 

x(x)  - 

2X-X* (x) 

'^l 

odd 

**(i) 

N, 

even 

*(x)  - 

2Z-z*(x) 

1 CM 

odd 

(42a) 


(42b) 
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i! 

li 

11 

c 


where  Is  the  total  number  of  reflections  of  the  characteristic 

at  X ■ 0 and  X,  and  N.  Is  the  number  of  reflections  at  z ■■  0.  We 
define 

“ 2X  ~ *2I~l+l/2*  ^ 1 ^ 21,  (43a) 

*J+l/2  " “ *2J-J+l/2’  < J < 2J.  (43b) 

With  Eq.  43,  the  definition  (Eq.  36)  of  R. , may  be  extended  to 
1<  1<  21,  1<  j < 2J. 

Points  yeD  %fhose  characteristics  lead  after  time  t to  a 
point  x*eD^  will  have 

x(y)e{x^,(x*,2Z-z*),  (2X-x*,z*),  (2X-x*,  2Z-z*)}.  (44) 

Equivalently,  let 

^Ij  'VJt  ~ (**(x)kRj^}  , 1^1^  I,  l^J^J,  l^k^I, 

1 < ^ < J,  (45a) 

and 

^lj;k4  ■ 1 < i<  1,  1 < j < J,  1 < k < 21, 

1<  / < 2J.  (45b) 

Then 

\j;ke  " ^lj;ke  ^ *^lj  ;2I+l-k,£  ^ *^lj  ;k,2J+l-^  ^ *^1J  ;2I+l-k,2J+l-£’ 
KKI,  1<J<J,  l<k<I,  1<^<J.  (46) 

Fluid  In  will  have  Its  x-veloclty  reversed  If  1 -t-  1 < k < 21 

and  Its  z-veloclty  reversed  If  J + 1 < ^ < 2J. 
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I 

I 


It  follows  from  Eq.  32  with  g ■ 0 and  the  assumption  that  u 
and  w are  constant  throughout  that  we  can  write 

I(l,j)  J(i.j) 

{x(x)|xcR  } - U U 

a-1  B-1  ^ 


where 


^lj;oB  = te'Cl.j.a),  ?'^(l,j,a)]  [n'd.j.B),  /(l.J.g)],  (47b) 


and  1(1, J)  Is  1 or  2,  according  to  whether  all  characteristics 


from  have  been  reflected  at  x 


0 and  X an  equal  or  unequal 


number  of  times,  respectively,  and  J(l,j)  Is  0,  1,  or  2,  according 
to  whether  all  characteristics  from  R^^j  have  left  the  computational 

grid,  or  whether  those  remaining  have  been  reflected  at  z - 0 an 
equal  or  unequal  number  of  times,  respectively.  To  find  1(1, j), 

J(i»j)»  5*.  h*,  we  do  the  following. 

We  first  construct 


J(1,J)  - 0; 

J(l,j)  - 1,  n"(l,j,l)  - z’j, 
- Z; 

z^j  * -Z  and  z^  A -Z:  J(1,J)  ■ 0; 

*Ij  * ^ il”(l,j,l)  - Z, 

and  n'’'(l,J,l)  - 2Z  + z^^; 


xj  jix/z  ij 


If 


z^j  ^ Z and  z”j  > Z: 


and  n (l,j,l) 


(48) 

(49a) 

(49b) 

(49c) 

(49d) 
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0 < < Z and  ^ 0;  J(l,j)  “ 1,  n ■ z^y 

and  n^(i,jtl)  “ ; 

0 < z^j  < Z and  z^j  < 0:  J(i,j)  “ 2,  n 

- 0,n'‘'(l,j,l)  - Zy.  n"(l,J,2)  - 2Z  + z"j,  and  n^(l 

-Z  < Zy  < 0 and  -Z  < z“j  < 0:  J(i,J)  - 1, 

- 2Z  + z~y  and  - 2Z  + z^^ . 

Next  we  find  the  unique  integer  m such  that 

*lj  “ ’^1-1/2  + “ijT  + 2m  X 
saClafles 

0 < < 2X 

and  the  unique  Integer  m^  such  that 

“ *1+1/2  “ij"  + 2“"’* 
satisfies 

0<  Xy  < 2X. 

If 

*ij  ^ *ij'  and 


(49e) 


J,2)  - 2Z; 
(49f) 


(A9g) 


(50a) 

(50b) 

(51a) 

(51b) 


*1J 
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s x"j:  1(1, j)  - 2,  rd.J.l)  - 

C"(i,j,2)  - x"j,  and  e‘*'(l,J,2)  - 2X.  (52b) 

Since  p,  u,  and  v are  assumed  consCanC  throughout  for 

each  1 and  j , we  have  for  ^he  mass  and  momentum  associated  with 
the  points  yeD^  such  that  x(y)E  Rj^,  1 < k < 21,  1 *:  £ S.  2J, 

I J 

■V-l  2 

1-1  j-1 
I J 

1-1  j-1 
I J 

*'zk£  "22  '"ij’ 

1-1  j-1 

where  |a|  Is  the  (Lebesgue)  measure  of  the  set  A.  From  Eq.  47, 

1(1. j)  J(l.j)  ^ 

^^ij;h^^  “ 2 2 l^ij;aB  ' 

a-1  S-1 

From  Eqs.  47b  and  36, 

'^lj;ae^Vl  - “ax|mln[5''(l,j,a),  x^^^/2]  ■ 

- max  [C"(i,j.oi).  *k_i/2]* 

X max  |min[n^(l,j,B), 

~ max  [n  (l.j.B),  2^_i/2]»  • ^5) 


(53a) 


(53b) 


(53c) 
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The  approximate  numerical  solution  of  the  hyperbolic  conservation 
law  is  completed  by  setting,  with  the  help  of  Eq.  46,  for  1 < 1 < 


I and 

1<  J < J, 

* 

* 

* * 

"ij- 

“21+1-1,  j ' 

“l,2J+l-J  “21+1-1, 2 J+l-J’ 

(56a) 

* 

* 

it  * 

^IJ " 

' *'xlj 

■ '^x2I+l-i,J 

+ \i,2J+l-j  ■ ''x2I+l-l,2J+l-j, 

(56b) 

* 

* 

* * 

' ‘'zlj 

‘^z2I+l-i,j 

■ ''zl,2J+l-J  " '^z2I+l-l,2J+l-j* 

(56c) 
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we  will  want  to  replace  Eq*  57b  by 


i 


where  u.  Is  the  value  of  u associated  with  cell  1.  When 

(Ax^) 

Is  sufficiently  small.  It  Is  most  convenient  to  approximate  S^(Aa)x^ 

by  a function  that  Is  also  constant  on  each  Interval  (x^  *ic+l/2^ 

and,  furthermore,  that  Is  zero  when  |k  - l|  > 1.  Thus,  we  shall 
write  approximately 

0 + 

*^1  *^1 

^±-1  Ax^  ^1  ^1+1 

where  c*  and  cj  are  to  be  found.  We  shall  determine  the  coeffi- 
cients by  requiring  that  the  first  three  moments  of  x be  equal  for 
the  expressions  on  the  left-  and  right-hand  sides  of  Eq.  59.  Re- 
ferring to  Eq.  57b,  we  easily  calculate 


Sx(Aa)Xi 


Ax 


11  Ax 


1-1 


I 

/ 


Sx(Aa)  X^dx  " AXj^, 


X S^(Aa)  Xj^dx  - x^  Ax^^  , 


J x^  S^(Aa)  Xj^dx  - (Xj^^  + -^(Ax^)^  + 2Aa)  Ax^, 


where 


(60a) 

(60b) 

(60c) 


X 1 

^ ■ 2^*l-l/2  *1+1/2^  • 


(60d) 


Equating  Eqs.  60a  through  60c  with  the  corresponding  moments  of 
the  right-hand  side  of  Eq.  59,  we  obtain  the  equations 
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? - 

Ik 


CiAXi 


0* 

CiAx^ 


+. 

CiAXi 


(63b) 


k - 1 - 1 


k - 1 


k - 1 + 1 


k - 1 > 1 


In  Eq.  20c,  we  observe  that  S^(Aa)  operates  on  functions 

that  have  been  continued  symmetrically  across  the  boundaries 
X - 0 and  X.  Thus,  If  one  were  to  represent  numerically  the  ef- 
fect of  S^(Aa)  at  cell  k on  a function  f that  has  the  value  f^ 

In  cell  1 and  Is  extended  to  a function  f according  to  the  pres- 
cription (Eqs.  20d  and  e),  the  result  would  be 


Ax.  Ik  1 Ax.  Ik  1' 


where 


2 ^1  «I-lor2  Sk  <;i-l 


Plk  - ( (^1  + cpAx^  1 - k - 1 


• (65) 


(Cj  + Cj)AXj  1 ■ k “ I 


All  the  quantities  In  Eq.  65  are  calculable  through  Eqs.  62  and 
63b,  with  the  convention  that 


^*0  “ ^*1/2  " 


^*1+1/2  " ^*1+1  ■ ^*1* 


(66a) 


(66b) 
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Similar  results  hold  for  the  approximation  of  the  operator 

S (Aa) . Let 
z 


*J  " 2^*J-1/2  *j+l/2^’ 


(67a) 


^^j+1/2  “ 2^^*j  ^*j+l^’ 


(67b) 


AZq  - “ ^*1* 


(67c) 


^*J+l/2  " ^*J+1  " ^*J- 


(67d) 


Then  the  effect,  evaluated  In  cell  t,  of  S (Aa)  operating  on  a 

z 

function  f that  has  the  value  f^  In  cell  J and  Is  extended  sym- 
metrically to  a function  f according  to  Eqs.  20d  and  20e,  Is, 

Act 

for =•  sufficiently  small. 


J 


where 


2 «j  «Jor2«tsj 


V 


(e^  + spAZj^  j - ^ - 1 
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0, 


+. 


and 


^ - J - 1 

I - J 

^ - j + 1 

1^  - j|  > 1 


(70) 


3da[dZj_i/2(d2j_i  + dzj  + dzj^j^)]”^. 

(71a) 

3da[dzj^l/2(d*j-l  + dzj  + AZj^^)]-\ 

(71b) 

-f 

1 - .j  - .j. 

(71c) 

When  da  Is  not  small  enough,  some  of  the  quantities  or 
Qjj  given  by  Eqs.  65  and  69  may  be  negative,  and  In  that  case  the 
approximate  representations  of  S (do)  or  S (da)  that  have  been 
given  above  %d.ll  lead  to  Instabilities  In  the  computation. 


Accordingly,  if,  for  example,  < 0,  It  will  be  necessary 
for  us  to  obtain  a better  representation  of  S^(da)x^  than  that  af- 
forded by  Eq.  59.  Referring  to  Eq.  57b,  we  look  at 


S,(4.)Xi 


x-x 


(4irda) 


172 


I 


1-1/2 


e-^ 


x-x 


1+1/2 


(72) 


when  X Is  In  the  kth  cell.  Consistent  with  the  accuracy  of  the 
numerical  quadrature  we  have  been  using,  we  replace  In  Eq.  63 
by  the  Integral  of  Eq.  72  over  the  kth  cell.  For  k > 1,  this  Is 
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? - 

Ik 


f 

lirAa  J 


>1/2  r"l-l/2  V/AAa,,  , 

e dC  dx 


»'4irAa  J J 

*k-l/2  *”*1+1/2 


p [ *k-l/2  ” *l+l/2\  _ p / *k+l/2  ~ *1+1/2 

2^  I I 2/Ea 


When  k < 1,  one  obtains 


(74) 


If  S^(Aa)  operates  on  a function  that  has  been  continued 

symmetrically  across  the  boundaries  x ■ 0 and  X,  the  result  may 
be  put  In  the  form  (Eq.  64)  where  now  P..  Is  given  approximately 
by 
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Similarly,  If  Q^j  calculated  by  Eqs.  69  through  71  Is  < 0, 
we  calculate  according  to  the  following  procedure.  Let 


*^+l/2  * ■ *j-l/2 


= /■  f 


/4Aa  , 
e dt  dz. 


(76) 


‘£-1/2 


'j+1/2 


When  Is  given  by  equations  similar  to  Eqs.  73  and  74. 

For  £ ■ J,  we  find 


- Az,  - 2 + 2pl 

J \2^ 


Az , 


(77) 


The  effect  of  S (Aa)  evaluated  at  cell  £ operating  on  a function 
z 

that  has  the  value  f ^ In  cell  J and  has  been  extended  symmetrically 

across  z * 0 according  to  Eqs.  20d  and  20e  may  be  put  In  the  form 
(Eq.  68),  where  Is  given  approximately  by 


Qj£-P 


1/2  ~ ^1+1/2  ) ( ^£+1/2  ~ ^l-H/2 ' 


- P 


l/Ea 

H ~ 

l/ET 


^£-1/2  ~ ^1-1/2 


+ pi 


2*^ 


+ P [ !ki£l_!izl/2  J p [ ^£+1/2  ^j-1/2 

l/Ea  I \ l/Ea 


. p,  :£iiZLlllt3Vi' 
2/Sa 


+ P 


^£-fl/2  ^1-t-l/2 

2/l^ 


1 < J < £ < J, 


(78a) 
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(78b) 


Q..  - A^.  - 2^  . 2P(-^V  P(M  - 2Pp 

J V U/aT/  I /a^T  / 


1 < j < j. 


(78c) 


In  the  numerical  solution  of  the  Stefan  problem,  we  follow 

the  prescription  laid  out  In  Eq.  20.  We  form  quantities  m"  , 
starting  with 


0 'V 

®lj  “ “ij* 


The  step  corresponding  to  Eq.  20c  In  which  Is  computed  Is  per- 

formed In  two  parts. 

First,  we  compute 


IJ  AXj^AZj 


and,  given  a small  positive  constant  we  check  whether  or  not 


(p“  - p ) < e . 

1<J<J  IJ  0 1 


The  constant  Is  what  terminates  the  solution  of  the  Stefan  prob- 
lem and  represents  an  allowable  margin  of  error  In  the  density  cal- 
culation. Suppose  Eq.  81  Is  first  satisfied  for  n - n..  Then  we 
set  ° 
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(Ay) 


xlj 


I 

AZj  2 


'kr 


n-O  k-1 


2 \ *1-1/2  ■*■  *1+1/2  " *k-l/2  ■ *k+l/2  ’ 


(87a) 


and 


(Ay)^lj  - A*l 


J 

2 I 

n-O  £^1 


,n+l/2j 


I 


i (*j-i 


1/2  ■*■  ^j+1/2  ■ ^Z-1/2  " ‘^+1/2, 


" “ Z/ 


(87b) 


Finally  we  determine  the  new  momenta  by  using  Eq.  24.  We 

calculate  quantities  y^  and  y*’  , starting  with 

*1J  ^ij 


and 


(Ay) 


zii’ 


max 

K i < I V.  <e, 
1<  j < J 


where  Is  a small  positive  constant,  we  set 


(88) 


(89) 
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'*XlJ  “ '^Xlj 

and 

'^zij  * 

Otherwise,  let 


max 

V - 1 ^ I V 

1 ^ J 


(90) 


(91) 


Mindful  of  the  sufficient  condition  for  stability  of  the  algorithm 
Eq.  24  as  given  In  Eq.  25,  we  set 


A - 1.1 


+ 

V 


> 


(92) 


Ay  ■ Aa/A, 


(93) 


and 


l.lv 


(94) 


Then  we  proceed  with  the  algorithm  Eq.  24  with  this  A and  Ay,  com- 
puting and  for  0 « p ^ p^  where  p^  Is  the  first  Integer 

for  which 


(Pq  + 1)  Ay  > y^. 


(95) 


and  y^  Is  a prescribed  positive  constant  with  e 
be  a small  number. 


considered  to 
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The  numerical  step  corresponding  to  Eq.  2Ac,  In  which 
and  found.  Is  done  In  two  parts.  First,  we  find 


jP+1/2  , y '^il 

xlj  Zj  Ax.  Az» 
£-1  ^ ^ 


(96a) 


and 


,P+l/2 

zlj 


^-1 


Az£ 


Then  we  compute 


(96b) 


I- 

ti 


,P+l/2 

xlj 


(1  - yP  + Y yP-^^/2 

^Ij  ^ ‘'xlj  ^ '^kl  ‘'xkj 


k-1 


and 


(97a) 


jP^^  “ (1  - v^^^)  yP  + Y P yP^^^^ 
zlj  ^ Ij  ^ Pzlj  ^ kl  Pzkj  • 

k-I 


(97b) 


Here  P^^^  and  are  found  In  a manner  reminiscent  of  the  deri- 

vation of  P and  Q above.  The  only  difference  Is  brought  about  by 

the  different  way  Is  obtained  from  In  Eq.  24  as  compared  with 

the  extension  of  O*'  to  s''  In  Eq.  20.  Specifically,  when  as 

calculated  by  Eqs.  62,  63,  and  65  Is  s 0, 
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and  when  Q. . calculated  by  Eqs.  69  through  71  Is  ^ 0 


If,  on  the  other  hand,  Eqs.  62,  63,  and  65  yield  a value  of  P , < 0 


la  given  by 


1+1/2 


(100a) 


(100b) 
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^ ^ "fD- 


X - X 


1+1/2 


X - X. 


X - X 


1-1/2 


1 a:  1 K I. 


(100c) 


And  If  Eqs.  69  through  71  give  < 0,  we  calculate  from  the 

alternative  formula 


- Q./  - 2P(  + 2P(  lbl/2  ^^+1/2 

V 2/ta  / \ 2/ta 


^ 2p  I ‘l-H/2  _ jp  / 

\ 2^0  / \ 2.^ 


1 S J K J,  1 A £ « J.  j^£, 


— 1-  _ p + 2P  -J- 

2/SoLr  ' \ i/ao  / \ /Aa  J 


. 1 K j K J. 


(101a) 


(101b) 


Finally I the  numerical  version  of  Eq.  23  Is 


■ I <i  - "I? 


(102) 
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4.  SAMPLE  CALCULATIONS 


We  shall  describe  some  calculations  that  we  have  made,  using  | 

the  numerical  representation  of  the  flow  described  In  Section  3.  A i 

listing  of  the  computer  program  appears  In  Appendix  A.  We  shall 

outline  the  results  of  the  calculations  one  by  one  and  examine  them  I 

for  any  Indications  they  offer  about  the  accuracy  of  the  code.  More  ] 

detailed  observations  about  the  deficiencies  of  the  present  program 
and  suggested  Improvements  will  be  given  In  Section  5. 

In  each  of  the  calculations,  we  have  found  It  worthwhile  to 
monitor  the  total  energy  as  a function  of  time.  In  general,  energy 
will  be  lost  In  the  part  of  the  algorithm  that  approximately  solves 
the  hyperbolic  conservation  laws  (Eq.  6).  This  follows  from  Eq.  31, 
which  expresses  the  fact  that  all  "collisions"  of  parcels  of  fluid 
are  Inelastic.  This  energy  loss  Is  not  Just  a function  of  the  time 
and  space  discretization  but  may  persist  even  In  the  continuous 
limit  and.  In  fact.  Is  Intimately  related  to  the  turbulent  or  non- 
turbulent  character  of  the  flow  (Ref.  2).  Nevertheless,  classical 
Invlscld  flows  conserve  energy,  and  by  examining  the  variation  of  * 

energy  with  the  time  we  may  get  an  Idea  either  of  the  degree  to 
which  Che  flow  la  turbulent  or  of  the  error  Involved  In  the  discrete 
approximation  of  the  flow.  The  energy  loss  due  to  the  discrete  ap- 
proximation essentially  has  two  sources:  a "diffusive"  energy  loss 
associated  with  the  finite  size  of  Ax  and  Az,  and  a "colllslonal" 
loss  associated  with  the  flnlteness  of  the  time  step  t and  the  pos- 
sibility of  different  characteristics  of  the  Boltzmann  equation 
(Eq.  28)  running  together  In  time  t. 

(The  algorithm  that  we  have  described  Is  not  guaranteed  to 
dissipate  energy,  since  energy  may  be  created  In  the  solution  of 
the  Stefan  problem  (Eq.  10)  and  the  elliptic  boundary  value  problem 
(Eq.  13).  However,  one  can  verify  that,  for  the  analytical  algo- 
rithm of  Section  2,  such  energy  production  does  not  exceed  the  col- 
llslonal loss  In  the  limit  as  t-M)  (Ref.  1).  Thus,  any  energy  In- 
crease that  takes  place  In  the  numerical  solution  of  the  discrete 
equations  must  reflect  the  discretization  error  Involved.) 

» 

EVOLUTION  OF  A LIOUID  INITIALLY  AT  REST  AND  WITH  A 
HIGHLY  DISTORTED  INITIAL  FREE  SURFACE 

Our  first  calculation  had  I - J - 10,  Ax^^  - AZj  - 1, 

T"0.1,  g*l,  and  pQ  ■ 1.  The  Initial  data  were 
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l<j<  2 

3<1<4,  1<J<10 
5<i<7,  l<j<5 
8<i<10,  l<j<3 
otherwise, 


(103) 


and  “0,  1<  1,  j < 10. 

culatlon  were  chosen  as  e > 10 

Aa  » 0.1,  as  defined  In  Eqs.  38,  67, 
The  calculation  was  run  until  tine  t 


Ocher  paraneters  for  the  cal- 
“ 10  Yq  “ 10,  and 

75,  81,  and  26,  respectively. 

- 14. 


Our  original  hope  In  choosing  the  Initial  liquid  domain 
as  given  In  Eq.  103  was  that  such  a highly  distorted  Initial  sur- 
face would  lead  to  wave  breaking  and  falling  over,  with  the  at- 
tendant formation  of  cavities.  However,  this  expectation  was  not 
borne  out  by  the  computational  output.  The  numerical  results  In- 
dicated a free  surface  for  which  the  vertical  coordinate  was  a 
single-valued  function  of  the  horizontal  coordinate.  In  this 
case  the  z-coordlnaCe  can  be  Identified  with  the  total  mass  In 
a column  of  fluid.  Table  1 gives  the  results  for  this  total  mass 
at  times  t ■ 0,  2,  4,  6,  and  8.  We  observe  that  the  free  surface 
appears  Co  oscillate  In  time,  with  the  oscillations  getting  pro- 
gressively smaller  as  time  Increases.  There  develops  a rather 
high  peak  of  the  free  surface  at  the  aide  walls  (1  1 and  10), 

a feature  that  has  dubious  authenticity.  Rather,  this  appears 
to  be  related  to  a defect  of  the  numerical  algorithm  of  Section  3 
with  regard  to  the  treatment  of  normal  components  of  velocities 
at  rigid  boundaries.  This  point  Is  discussed  more  fully  In  Sec- 
tion 5,  where  we  also  suggest  an  Improvement. 

In  Table  2,  we  show  Che  kinetic  energy  (K£),  potential 
energy  (PE) , and  total  energy  (E)  of  the  fluid  system  as  a func- 
tion of  time.  We  note  that  there  Is  an  Initial  Increase  In  en- 
ergy. Of  course,  this  Is  a spurious  effect  and  Indicates  the 
size  of  the  discretization  error  made.  As  time  progresses,  the 
energy  appears  to  decay.  We  have  pointed  out  above  that,  for  a 
nonturbulent  flow,  such  decay  Is  not  a property  of  the  exact 
solution  but  reflects  the  error  In  the  finite  representation  of 


- 44  - 


i 

I 


t! 

fl 

{1 

IJ 


I 


Q 


THE  JOHNS  HOmiNS  UNIVENSITV 

APPUED  PHYSICS  LABORATORY 
LAUREL  Maryland 


r 


I 

I 


Table  1 

Total  mass  in  column  i as  a function  of  t. 
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Table  2 

Kinetic  energy,  potential  energy,  and  total  energy  as  functions  of 
t,  for  a highly  distorted  initial  surface. 
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the  evolutionary  equations.  In  this  problem,  because  of  the  ab- 
sence of  any  observed  breaking  or  development  of  other  pathologies 
In  the  flow,  we  are  inclined  to  discount  the  presence  of  turbu- 
lence. The  energy  decay  Is  such  that  we  would  tend  not  to  give 
much  credence  to  the  quantitative  computer  results  after  time  > 8, 
even  apart  from  the  unpersuasive  character  of  the  data  at  the  side 
walls  for  earlier  times. 

COMPARISON  WITH  LINEAR  THEORY 

The  algorithm  presented  in  this  report  Is  decidedly  ineffi- 
cient when  It  comes  to  solving  linear  wave  problems.  First,  in 
the  linear  regime  the  problem  ceases  to  be  a free-boundary  problem 
in  any  Important  respects,  and  methods  based  on  Green's  function 
for  the  unperturbed  domain  are  more  effective.  Also,  a special 
burden  Is  placed  on  the  size  of  the  numerical  mesh,  as  It  must  be 
fine  enough  to  resolve  the  linear  displacements  of  the  free  sur- 
face, and  yet  the  assumption  of  linearity  requires  that  this  be 
only  a small  portion  of  the  vertical  extent  of  the  computational 
domain.  Nevertheless,  It  is  Incumbent  on  us  to  compare  the  results 
of  a calculation  based  on  our  algorithm  with  a known  solution,  and 
in  this  respect  a linear  problem  naturally  comes  to  mind. 

The  linear  solution  we  compare  with  Is  the  wave  whose  sur- 
face height  Is  given  by 


z(x,t)  ■ 5 + cos  ut  cos  0 < X iC  20,  t > 0,  (104a) 


where 


2 

u 


gk  tanh  kh,  g 


4,  h 


(104b) 


Our  computation  took  place  on  the  mesh 


Ax^  - 1,  1«  1 « 20, 


(105a) 


0.5 


1 « j * 3 
4 < J < 11 


(105b) 
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Initial  values  of  corresponding  to  the  Initial  profile 


2(x,0)  - 5 + cos  , 


(106a) 


were  provided,  and  the  Initial  values  of  the  momenta  were 


‘'xlj  “ *'zlj  “ 0.  1 ^ i ^20*  1 ^ J ^ 


(106b) 


These  momenta  were  consistent  with  Che  fact  Chat  the  fluid  whose 
free  surface  Is  given  by  Eq.  104a  Is  at  rest  at  t > 0.  Other 

relevant  constants  for  the  calculation  were  c ■ 10  ■ 0.005, 

_3 

■ 10  » Yq  " ^0  " program 

was  run  up  to  time  t ••  5.9. 

As  we  observed  In  our  discussion  of  the  first  example,  the 
z-coordlnate  of  the  free  surface  can  be  Identified  with  the  total 
mass  In  a column  of  fluid.  Table  3 shows  values  of  the  total  mass 

In  each  column  for  times  t *■  0,  0.5,  1.0,  1.5,  2.0,  and  2.5.  These 

are  to  be  compared  with  the  values  (Eq.  104)  predicted  by  linear 
theory  and  given  In  Table  4.  As  In  the  first  example,  we  appear 
to  get  an  accumulation  of  fluid  at  the  side  walls.  This  accumu- 
lation becomes  noticeable  after  t ■ 1.5.  By  t - 1.0,  we  appear 
to  be  departing  from  the  monotone  dependence  on  x predicted  by 
linear  theory.  This  departure  seems  to  commence  at  the  side  walls. 
The  data  for  times  t > 2.5  Indicate  greater  departures  from  the 
linear  solution.  It  Is  possible  that  genuine  nonlinear  effects 
should  arise,  since  the  Initial  height  varies  from  6 to  4 as  x 
varies  from  0 to  20.  If  the  program  we  have  described  were  con- 
sidered to  be  a final  product,  we  would  be  well  advised  to  con- 
sider this  point. 

Table  5 gives  the  kinetic  energy,  potential  energy,  and 
total  energy  of  the  fluid  as  a function  of  time.  A peculiar  fea- 
ture is  that  Initially  the  potential  energy  varies  only  slightly, 
undergoing  a slow  steady  decay.  Of  course,  this  behavior  Is  not 
consistent  with  the  linear  theory.  It  Is  only  after  about  t - 4.4 
that  the  total  energy  remains  essentially  constant.  The  kinetic 
energy  oscillates  from  one  time  to  the  next  and  also  undergoes  a 
slower  oscillation,  which  has  a peak  with  center  around  t * 2.5 
and  a low  with  center  around  t ~ 4.8.  Analysis  of  the  mass  totals 
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Table  3 

Total  mass  in  column  i as  a function  of  t. 


1 

0 

0.5 

5.997 

5.956 

5.972 

5.913 

5.924 

5.861 

5.853 

5.795 

5.760 

5.708 

5.649 

5.609 

5.523 

5.496 

5.383 

5.356 

5.233 

5.208 

5.078 

5.073 

4.922 

4.932 

4.767 

4.774 

4.617 

4.627 

4.478 

4.515 

4.351 

4.396 

4.240 

4.274 

4.147 

4.191 

4.076 

4.124 

4.028 

4.077 

4.003 

4.062 

5.875 

5.748 

5.715 

5.651 

5.581 

5.504 

5.402 

5.290 

5.158 

5.056 

4.951 

4.811 

4.680 

4.556 

4.514 

4.424 

4.326 

4.248 

4.197 

4.207 


5.759 

5.492 

5.552 

5.437 

5.401 

5.328 

5.263 

5.186 

5.091 

5.032 

4.958 

4.861 

4.777 

4.647 

4.592 

4.564 

4.599 

4.483 

4.379 

4.442 


5.614 

5.184 

5.338 

5.191 

5.197 

5.135 

5.108 

5.057 

5.013 

4.986 

4.945 

4.899 

4.873 

4.766 

4.692 

4.648 

4.735 

4.876 

4.688 

4.847 


5.406 

4.889 

5.123 

4.927 

5.021 

4.930 

4.948 

4.909 

4.911 

4.914 

4.909 

4.905 

4.909 

4.886 

4.816 

4.757 

4.819 

4.999 

5.217 

5.554 
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Table  4 

Position  of  the  free  surface  for  a linear  wave  at  the  center  of  cell 
i as  a function  of  t. 
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Table  5 

Kinetic  energy,  potential  energy,  and  total  energy  as  functions  of 
t,  for  a slightly  distorted  initial  surface. 


59.82 
52.56 

54.82 
53.16 

51.54 

52.88 

49.89 

50.70 

48.24 
46.46 

45.52 

43.50 
41.41 
42.77 
40.20 

41.37 
40.27 
40.26 
39.95 

40.55 
40.01 

41.24 

40.50 

41.38 
41.10 

41.71 
41.74 
42.73 

43.53 

44.71 


t 

— 

KE 

1 

PE 

E 

0 

1021 

1021 

0.1 

13.78 

1018 

1032 

0.2 

31.08 

1018 

1049 

0.3 

28.35 

1017 

1045 

0.4 

44.76 

1016 

1061 

0.5 

31.58 

1016 

1047 

0.6 

43.67 

1015 

1058 

0.7 

37.55 

1014 

1051 

0.8 

45.65 

1013 

1058 

0.9 

42.59 

1012 

1055 

1.0 

46.61 

1011 

1058 

1.1 

47.68 

1010 

1058 

1.2 

49.37 

1010 

1059 

1.3 

51.94 

1009 

1061 

1.4 

52.13 

1009 

1061 

1.5 

54.05 

1008 

1062 

1.6 

54.23 

1008 

1062 

1.7 

55.65 

1008 

1063 

1.8 

56.39 

1008 

1064 

1.9 

56.97 

1007 

1064 

2.0 

57.87 

1007 

1065 

2.1 

58.27 

1007 

1065 

2.2 

59.22 

1006 

1065 

2.3 

59.46 

1005 

1065 

2.4 

60.61 

1005 

1065 

2.5 

60.97 

1004 

1065 

2.6 

60.79 

1003 

1064 

2.7 

61.86 

1003 

1065 

2.8 

58.82 

1004 

1063 

2.9 

54.60 

1004 

1058 
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for  the  different  columns  of  fluid  and  different  times  Indicates  a 
surface  that  Is  relatively  flat  at  t » 2.5.  By  way  of  comparison, 

the  linear  theory  predicts  a surface  that  Is  flat  at  t = . 

At  this  time  the  kinetic  energy  would  be  a maximum,  and  It  would 

then  decrease  to  0 at  t “ — . In  our  calculation,  the  "minimun" 

0) 

kinetic  energy  Is  about  two-thirds  the  "maximum"  value.  Thus,  from 
the  point  of  view  of  location  of  the  free  surface  and  period  of  os- 
cillation, our  calculation  gives  results  as  good  as  can  be  expected 
for  the  grid  we  have  used;  but  from  the  point  of  view  of  energy 
balance,  the  picture  Is  not  as  satisfactory. 

Another  feature  of  the  flow  that  can  be  compared  with  the 
prediction  of  linear  theory  Is  the  "pressure."  We  do  not  compute 
a pressure,  but  In  the  Interior  of  the  region  of  flow  the  quantity 
2 

— j v takes  the  place  of  the  pressure  for  classical  hydrodynamic 

T 

flows  (Ref.  1).  In  the  linear  theory,  the  pressure  should  be  es- 
sentially the  hydrostatic  pressure,  or  p.g  times  the  distance  be- 

14 

low  the  free  surface.  In  Table  6 we  give  Vj^j^  and  ^ m^^^^  j as 

j“l 

a function  of  time.  After  some  Initial  oscillation,  the  values 
of  V settle  down  around  t ° 1.  Thereafter,  they  appear  to  agree 
very  well  with  the  values  predicted  by  the  linear  theory. 

(Note  that  Eqs.  12  and  lOd  and  e Imply  that,  for  the  ana- 
lytical algorithm,  Vvn  ■ 0 for  xe9D.  This  is  not  true  of  the 

2 

the  pressure,  and.  In  fact,  — y v will  differ  from  the  pressure  In 

a layer  of  thickness  -s—  g near  z * 0 where  v will  jump  frOm  0 to 

z 

“ "2“  8 

As  we  observe  from  Eq.  26,  for  the  computational  scheme 
described  In  this  report,  v will  be  a rough  measure  of  the  compu- 
tational time  taken  In  solving  the  Stefan  problem  (Eq.  10).  This 
Is  the  longest  part  of  the  calculation,  and  accordingly  we  may  ex- 
pect the  computational  time  overall  to  Increase  with  the  size  of 
the  expected  values  of  v,  when  the  method  of  solution  Is  the  one 
we  have  used  to  date.  In  the  next  section,  we  will  discuss  more 
efficient  ways  of  computing  v. 


1 
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Table  6 

Total  mass  between  x = 10  and  x = 1 1 and  vi  as  functions  of  t. 


t 

14 

j-1 

^11.1 

1 

14 

^ “ll  i 
J-1 

"ll.l 

■ 

14 

^ “ll  J 
j-1 

^1.1 

0.1 

4.919 

0.12901 

2.1 

4.940 

0.09494 

4.1 

0.09830 

0.2 

4.921 

0.10557 

2.2 

4.933 

0.09601 

4.2 

0.09385 

0.3 

4.923 

0.07182 

2.3 

4.926 

0.09381 

4.3 

mm 

0.10251 

0.4 

4.928 

0.13092 

2.4 

4.918 

0.09579 

4.4 

m 

0.09087 

0.5 

4.932 

0.05875 

2.5 

4.909 

0.09476 

4.5 

4.841 

0.10141 

0.6 

4.935 

0.12553 

2.6 

4.900 

0.09416 

4.6 

4.885 

0.09525 

0.7 

4.940 

0.07098 

2.7 

4.890 

0.10167 

4.7 

4.940 

0.09883 

0.8 

4.944 

0.12000 

2.8 

4.880 

0.09101 

4.8 

4.993 

0.09677 

0.9 

4.948 

0.07786 

2.9 

4.871 

0.08565 

4.9 

5.031 

0.10034 

1.0 

4.951 

0.10537 

3.0 

4.861 

0.11353 

5.0 

5.051 

0.09633 

1.1 

4.953 

0.09245 

3.1 

4.850 

0.07945 

5.1 

5.064 

0.10168 

1.2 

4.955 

0.09562 

3.2 

4.840 

0.10154 

5.2 

5.074 

0.09318 

1.3 

4.957 

0.09915 

3.3 

4.829 

0.09005 

5.3 

5.083 

0.09861 

1.4 

4.958 

0.09356 

3.4 

4.818 

0.09396 

5.4 

5.092 

0.09601 

1.5 

4.958 

0.09809 

3.5 

4.807 

0.09789 

5.5 

5.103 

0.09886 

1.6 

4.957 

0.09404 

3.6 

4.797 

0.08945 

5.6 

5.117 

0.09541 

1.7 

4.956 

0.09700 

3.7 

4.789 

0.10097 

5.7 

5.133 

0.09870 

1.8 

4.954 

0.09573 

3.8 

4.783 

0.09340 

5.8 

5.148 

0.09665 

1.9 

4.950 

0.09536 

3.9 

4.779 

0.09755 

5.9 

5.164 

0.09752 

2.0 

4.945 

0.09534 

4.0 

4.778 

0.09741 

•TBr- 
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Since  the  linear  flow  is  Irrocatlonal,  a further  test  of 
the  accuracy  of  our  algorithm  would  be  a check  on  the  vorticlty 
of  the  computed  flow.  The  same  sort  of  test  might  also  be  per- 
formed on  the  data  of  the  first  example,  since  we  would  expect 
that  flow  to  be  Irrotatlonal  in  the  absence  of  breaking.  We  have 
not  examined  the  vorticlty  for  these  flows  in  detail  because  the 
preliminary  nature  of  our  results  does  not  seem  to  warrant  It  at 
this  time. 

COLLISION  OF  STREAMS  WITH  JET  FORMATION 

We  performed  three  runs,  with  different  computational  meshes 
and  time  steps,  for  the  flow  corresponding  to  the  initial  conditions 


p(x,z) 


0<x<2,  0 < z < 7 


0<x<2,  z>7 
2<x<5,  0 < z < 2 


2 < X < 5,  z > 2 


(107a) 


(pu) (x,z) 


-10p(x,z) 


0 < z < 2 


z > 2 


(107b) 


-10p(x,z) 


0<x<2,  z>2 


(pw) (x,z) 


0<x<2,  0<z<2.  (107c) 


X > 2 


-5  '2 

For  all  three  runs,  we  had  g ■ 1,  p.  1,  e - 10  , e,  - 10  , 

-3  u 1 

^ fl  n j __  ..  « 


Ej  ■ 10  , and  Yq  ■ 10.  Otherwise,  we  had  for  run  1 


- 54  - 
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T ■ 0.1,  “ 1,  AZj  "1,  I - 5,  J ■ 10,  Aa  ■ 0.1,  (108a) 

and  we  ran  Che  problem  10  time  steps;  for  run  2 

T - 0.05,  Ax^  - 0.5,  AZj  - 0.5,  I - 10,  J - 20,  Ao  » 0.05,  (108b) 

and  the  problem  was  run  20  time  steps;  for  run  3 


T - 0.025,  Ax^  - 0.25,  AZj  - 0.25,  I - 20,  J - AO,  Aa  - 0.025, 

(108c) 

and  we  computed  the  flow  for  40  time  steps. 

Tables  7,  8,  and  9 record  the  kinetic  energy,  potential  energy, 
and  total  energy  as  functions  of  time  for  the  three  runs.  Gener- 
ally we  observe  that  the  energy  tends  to  be  higher  at  a given  time 
for  the  run  with  the  finer  computational  mesh.  This  Is  In  accor- 
dance with  our  observation  at  the  beginning  of  this  section  that 
the  finite  grid  leads  to  a spurious  energy  loss  through  diffusion 
and  collisions.  However,  beyond  this  energy  loss  there  appears  to 
be  an  energy  loss  for  t < 0.3  that  Is  not  related  Co  the  finite  grid 
spacing  but  that  may  reflect  the  presence  of  turbulence  In  the  flow. 
After  about  t ■ 0.3  the  slower  diminution  of  energy  observed  may  be 
due  primarily  to  the  error  Inherent  In  the  discretization.  (However, 
the  exact  solution  would  still  be  expected  to  exhibit  energy  loss 
associated  with  the  collapse  of  cavities  after  t - 0.3.)  As  we  ob- 
served In  Che  second  example.  In  all  three  runs  the  potential  energy 
varies  slowly. 

The  three  runs  were  compared  for  their  consistency  In  depict- 
ing the  free  surface  at  a given  time.  The  hope  Is  that  one  can  get 
a measure  of  the  error  In  a computation  by  examining  the  dependence 
of  the  output  on  the  mesh  size.  Of  course,  agreement  of  calculations 
and  the  demonstration  of  their  convergence  says  nothing  about  what 
they  converge  to.  That  Is  a task  for  the  theory.  We  chose  to  be 
rather  crude  In  plotting  the  free  surfaces  obtained  In  order  not  to 
give  the  numerical  results  any  particular  advantage.  Our  criterion 
for  drawing  a free  surface  Is  as  follows:  If  p..  as  computed  In 
1 

Eq.  39  la  > -r,  the  cell  Is  Included  In  the  "water"  region;  If 
1 ^ 

< -J,  the  cell  Is  In  the  "vacuum"  region. 
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Table  7 

Kinetic  energy  as  a function  of  t for  three  runt. 


0 

.025 
.050 
.075 
.100 
.125 
.150 
.175 
.200 
.225 
0.250 
0.275 
0.300 
0.325 
0.350 
0.375 
0.400 
0.425 
0.450 
0.475 
0.500 


Run  1 Run  2 Run  3 


1000  1000 


483.74 


361.23  404.73 


349.13 


274.86  254.41 


218.99 


129.25  170.10 


144.44 


92.54  130.18 


122.11 


74.51  114.61 


1000 

628.12 

575.56 

519.95 
481.68 
431.86 

392.15 
326.13 

301.95 
262.88 
238.81 
221.39 
209.67 
200.55 
193.09 
186.75 
180.31 
174.73 
171.17 

168.95 

167.16 


.525 

.550 

0.575 

.600 

.625 

0.650 

0.675 

0.700 

0.725 

0.750 

0.775 

.800 

.825 

.850 

0.875 

.900 

.925 

.950 

0.973 

1.000 


: 1 

Run  1 

Run  2 

Run  3 

109.11 

164.04 

158.33 

65.99 

103.28 

149.79 

138.34 

96.00 

130.86 

116.01 

60.91 

87.38 

107.38 

101.45 

76.43 

95.54 

87.07 

52.89 

59.65 

73.58 

64.65 

53.32 

59.69 

55.29 

44.53 

48.28 

51.83 

47.85 

42.93 

44.25 

43.40 

39.36 

39.81 

42.29 

42.47 
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0 

.025 

.050 

.075 

.100 

.125 

.150 

.175 

.200 

0.225 

0.250 

0.275 

.300 

.325 

.350 

0.375 

0.400 

0.425 

0.450 

0.475 

0.500 


Table  8 

Potential  energy  as  a function  of  t for  three  runs. 


Run  1 

Run  2 

Run  3 

55 

55 

55 

53.70 

52.78 

53.51 

53.38 

51.36 

52.70 

53.25 

53.11 

52.66 

52.97 

52.84 

51.89 

52.59 

52.69 

52.52 

52.48 

52.34 

52.18 

52.48 

52.32 

52.01 

51.82 

52.13 

51.64 

51.45 

53.05 

51.89 

51.21 

50.96 

51.59 

50.74 

50.51 

53.51 

51.20 

50.24 

0.525 

0.550 

0.575 

0.600 

0.625 

0.650 

0.675 

0.700 

0.725 

0.750 

0.775 

0.800 

0.825 

0.850 

0.875 

0.900 

0.925 

0.950 

0.975 

1.000 


Run  1 


53.80 


53.95 


52.80 


52.12 


51.62 


Run  2 


50.73 


50.06 


49.25 


48.68 


48.82 


47.81 


47.53 


47.37 


47.17 


46.88 


Run  3 


48.33 

47.39 

46.66 
46.19 
45.80 
45.46 

45.12 
44.75 

44.29 
44.21 
44.05 
43.84 

43.66 

43.40 

43.12 

43.30 
43.40 
43.52 
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Run  1 


55 


326.75 


181.72 


145.54 


Table  9 

Total  energy  as  a function  of  t for  three  runs. 


536.52 


412.59  I 457.43 


401.79 


128.02 


307.00 


271.46 


222.42 


196.56 


182.07 


173.70 


165.81 


Run 

3 

1055 

681, 

,82 

629. 

,07 

519. 

,95 

534, 

,93 

484, 

,97 

445. 

,12 

378. 

,97 

354. 

,64 

315. 

,40 

291. 

.16 

273. 

,57 

261. 

.68 

252. 

,37 

244. 

,72 

238. 

,20 

231. 

,52 

225. 

,69 

221. 

,92 

219. 

.46 

217. 

,40 

nili] 


Run  1 Run  2 


159.84 


119.79 


114.86 


105.69 


96.65 


90.98 


153.34 


145.25 


136.07 


125.25 


107.46 


100.85 


95.65 


90.10 


86.69 


Run  3 


213.86 

207.51 
198.13 
185.73 

177.52 
162.20 
153.19 
146.91 
140.66 
131.81 

117.87 
108.85 
103.75 

99.12 

95.50 

91.25 

87.37 

86.71 

85.69 

85.99 
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Figure  1 shows  the  Initial  water  surface.  Figures  2 through 
31  depict  computed  water  surfaces  for  various  times  and  various 
runs.  In  general,  the  agreement  among  the  figures,  especially 
those  for  runs  2 and  3,  Is  rather  good.  Although  we  cannot  have 
too  much  faith  In  the  runs  for  later  times  because  of  the  loss  of 
energy.  It  Is  still  not  unreasonable  to  expect  them  to  exhibit 
correctly  some  of  the  qualitative  features  of  the  flow.  Thus,  we 
may  expect  that  for  the  actual  flow  a cavity  appears  In  the  left 
Interior  around  t - 0.1,  and  that  by  t - 0.2  a jet  has  struck  the 
right  wall  X ••  5 and  a cavity  has  been  formed  there.  The  Interior 
cavity  at  the  left  disappears  between  t ■ 0.5  and  t ■ 0.8,  and  the 
cavity  at  the  right  closes  In  at  the  wall  around  t ■ 0.7.  By  t » 
0.8  the  larger  cavities  have  closed  In.  One's  Intuition  might 
lead  one  to  expect  the  jet  to  bounce  off  the  right  wall  to  create 
a leftward  moving  jet.  Indeed,  we  see  a hint  of  such  behavior  In 
Figs.  10  and  13.  It  Is  possible  that  the  program  has  suppressed 
this  tendency  by  allowing  fluid  to  accumulate  at  the  right-hand 
wall  Instead,  In  a manner  reminiscent  of  the  computed  flow  for  the 
first  two  examples  above. 
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5.  LIMITATIONS  AND  IMPROVEMENTS  OF  THE 
PROGRAM 


One  of  the  most  serious  limitations  of  the  program  Is  In 
the  solution  of  the  hyperbolic  conservation  laws  (Eqs.  6 through 
9).  We  observe  from  Eqs.  39  and  42  through  56  that  our  algorithm 
allows  a diffusion  of  mass  across  the  computational  grid  In  addi- 
tion to  the  convective  flow  described  by  the  conservation  laws, 
and  this  reduces  the  clarity  of  the  delineation  of  the  free  sur- 
face. There  are  relatively  simple  remedies  available  at  slight 
cost  In  computational  effort.  One  such  remedy  Is  along  the  follow- 
ing lines:  during  the  convective  process,  when  we  move  parcels  of 
fluid  from  one  cell  to  another,  we  may  associate  with  the  fluid 
so  moved  not  only  a mass  and  momentum,  but  also  a center  of  mass. 
Thus  at  each  time  we  may  assign  a center  of  mass  to  the  fluid  In 
each  cell.  And  we  may  consider  the  fluid  In  the  cell  1,J,  with 
center  of  mass  *oij^  reside  In  a rectangle  of  area 


A - 4 mln(xQ^j  - ’'l+l/2  " *01J^  ' 'j-1/2' 

*J+l/2  ■ *01J^' 


unless  the  ratio  of  the  mass  In  the  cell,  j i A exceeds  p^. 

In  that  case  we  may  consider  m. . to  be  uniformly  distributed  over 

“l1  ^ 

cell  l,j  with  density  , as  we  have  done  heretofore.  By  such 

a procedure,  we  can  limit  the  diffusion  of  mass  due  to  the  finite 
cell  size.  As  a practical  matter,  we  have  found  this  spurious 
diffusion  of  mass  to  be  greatest  In  the  case  where  g ••  0.  When 
g > 0 we  have  observed,  not  unexpectedly,  that” the  gravity  tends  ” 
to  stabilize  the  free  surface,  which  Is  usually  confined  to  one 
or  two  computational  cells  In  thickness. 


The  novelty  of  our  approach  to  hydrodynamics  lies  In  the 
replacement  of  the  usual  divergence  condition  on  the  velocity  by 
the  constraint  p ^ p^.  That  part  that  deals  with  the  hyperbolic 

conservation  laws  Is  not  new,  at  least  from  a computational  point 
of  view.  It  may  be  that  other  numerical  work  on  such  conservation 
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laws  Is  more  satisfactory  than  our  own  treatment,  and  that  problems 
such  as  the  mass  diffusion  just  referred  to  have  already  been  ade- 
quately handled  in  other  Investigations.  Work  currently  In  pro- 
gress by  M.  Y.  Hussalnl  (Ref.  5)  uses  our  treatment  of  the  density 
constraint  In  a three-dimensional  Incompressible  flow  and  solves 
the  hyperbolic  conservation  laws  using  a MacCormack  "higher-order" 
hyperbolic  solver  (Ref.  6). 

In  the  examples  reported  In  the  last  section,  we  noted  the 
need  for  an  Improved  treatment  of  the  velocities  at  the  rigid 
boundaries.  Characteristically,  we  find  a rather  large  outward 
normal  velocity  at  the  cells  adjacent  to  the  rigid  boundary.  This 
sort  of  behavior  Is  encouraged  by  our  numerical  representation  of 
Eq.  27  as  Eq.  87.  We  may  expect  a more  satisfactory  treatment  by 
regarding  the  right-hand  side  of  Eq.  27  as  an  Integral  over  all 

x'eR^  and  6*'(x')  extended  s)nmnetrlcally  across  the  rigid  boundary 
dD. 


The  determination  of  v can  Itself  be  made  more  efficient 
than  the  method  used  In  Section  2,  where  a Stefan  problem  was 
solved  until  steady  state  was  reached.  For  example,  one  may  make 
use  of  the  monotone  dependence  of  the  solution  of  the  steady-state 

*\j 

Stefan  problem  on  p,  and  also  of  the  fact  that  this  solution  may 
be  obtained  by  solving  a succession  of  N steady-state  Stefan  prob- 

N ^ 

lems  with  Initial  data  0,  ^ p^  = p (Ref.  7),  to  obtain 

1“1 

directly  a lower  approximation  to  v,  with  the  remainder  of  v be- 
ing determined  Iteratively.  For  example,  this  would  be  desirable 
In  the  solution  of  problems  In  water  of  great  depth,  where  v,  be- 
ing proportional  to  the  pressure,  would  get  quite  large. 


P.ef.  S.  H.  Y.  Hussalnl  (private  communication). 

Ref.  6.  R.  W.  MacCormack,  "An  Efficient  Numerical  Method 

for  Solving  the  Time-Dependent  Compressible  Navler-Stokes  Equations 
at  High  Reynolds  Number,"  Comput.  Appl.  Math.,  Vol.  18,  1976,  p.  49. 

Ref.  7.  J.  C.  W.  Rogers,  "Steady  State  of  a Nonlinear  Evo- 
lutionary Equation,  Semlnalres  IRIA,  Analyse  et  Controle  de  Systemes, 
1978. 
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An  example  of  small  Improvements  that  might  be  made  In  the 
program  Is  the  following:  at  present,  we  only  solve  the  constraint 
(Eq.  1)  approximately,  getting  p * where  Is  given  In 

Eq.  67.  Thus  we  would  expect  the  density  computed  In  the  liquid 
domain  at  each  time  to  exceed  p^  by  a small  amount  proportional 

to  e^,  and  this  In  turn  should  lead  to  some  "settling"  of  the 

liquid  (In  the  direction  of  the  gravitational  force) . This  situ- 
ation can  be  ameliorated  by  solving  the  Stefan  problem  with  p^  re- 

placed  by  Pq 2 the  definition  of  the  function  f in  Eq.  10b, 

and  replacing  the  test  (Eq.  67)  by  the  test 


max 

1 ‘ j 


For  the  future,  the  first  thing  we  would  like  to  do  Is  to 
improve  the  treatment  of  velocities  at  the  rigid  boundary,  espe- 
cially the  numerical  representation  of  Eq.  27.  Beyond  that,  we 
are  thinking  of  making  the  code  applicable  to  the  computation  of 
Internal  waves  In  stratified  fluids.  This  would  require  only  a 
relatively  modest  addition  to  the  program  as  it  now  stands  (Ref.  1) . 
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Appendix  A 

PROGRAM  DESCRIPTION  AND  LISTING 

The  following  water  wave  program  was  written  for  the  opti- 
mizer and  checkout  PL/I  compilers  and  executed  on  ar  IBM  360/91 
computer  at  the  Frank  T.  McClure  Computing  Center  of  APL  (Ref.  8). 

Originally,  the  program  was  written  as  one  long  program, 
but  we  found  that  initial  conditions  were  easier  to  program  in- 
line, rather  than  read  in  as  input  data,  so  the  program  was  broken 
into  various  sections. 

The  main  procedure  first  states  various  constants  for  a 
given  run.  It  also  tests  certain  conditions  for  convergence,  when 
to  stop  and  when  to  print  answers,  and  when  to  write  on  a disk  in 
order  to  restart  or  continue  a problem  at  a future  time. 

Procedure  INITAL  computes  the  x-  and  z-coordinates  of 
points  in  the  extended  computational  grid. 

Procedure  PSAQS  computes  the  matrix  elements  that  simulate 
the  effect  of  diffusion  in  the  x-  and  z-dlrections  by  transforming 
quantities  of  mass  and  momentum  in  each  computational  cell  into 
new  values  through  multiplication  by  the  appropriate  matrices. 

Procedure  MASMON  computes  the  effect  of  convection  for  a 
time  step  on  the  values  of  mass  and  momentum  in  each  computational 
cell. 


Procedure  DENSTY  computes  a new  set  of  masses  for  each  com- 
putatiohal  ceir“at~tlie~en3  of  each  time  step  by  satistying  the 
constraint  on  the  density. 

Procedure  MOMEN  computes  the  final  anounts  of  momentum  in 
each  computational  cell  at  the  end  of  each  time  step. 

Procedure  PRNTAL  does  what  its  name  implies  — it  prints  out 
the  desired  information. 


Ref.  8.  "The  Frank  T,  McClure  Computing  Center  User's 
Guide,"  APL/JHU  BCS-1-92,  1 Sep  1978. 
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EXCEPT  MMERE  NOTEO.  all  VANUSlES  AKE  single  P«£CI:»10iy 

float  binary  numbers. 

THE  exception  to  The  rule  are  Variables  i.imaxeimi.imz...,, 
nHICH  follow  the  normal  naming  CUNVENTIUNS. 

IF  A VARIABLE  IS  OIMfNSIONEUt  IT  WILL  Bt  NUTfD  BY  PaRENThLSES, 


variables 

where 

USED 

OtSCRIRTION 

DA 

main.  PSAOSt 

densty.  momen 

STEP  SIZE  OF  InOEPENOEmT  VARIABLE 
OF  EQUATION  10  A, 

LABELED  "Aa"  IN  THE  TtXT. 

DEBUG 

Main,  densty 

MOMEN 

BITiniA  TRUTM/FALSE  SWITCH  TO 
pHooucE  debug  UUT^UT, 

D6 

MOMEN 

STEP  SIZE  OF  InOEPENOEmT 
variable  y OF  EQUATION  22  A. 

DT 

main.  maSmOn. 
OENSTy.  momen 

TIME  STEP,  labeled  aS  •' t •• 

IN  THE  text. 

DX (20) 

main.  PSAOS. 

INITAL.  MASMONi 
OENSTy.  momen 

widths  OF  CELLS  Rii 

IN  EQUATION  36. 

DZ (AO) 

main.  PSAOS. 

INITal.  m*smon. 
OENSTy.  momen 

HEIGHTS  OF  CELLS  Rjj 

IN  EQUATION  36. 

EPS 

main,  PPNTALt 
MASMON 

small  cut-off  to  KEEP  FROM 

dividing  by  zero  in  equation  38. 

EPSl 

MAIN,  densty 

small  PakamETER  which  determines 
when  DENSITY  Constraint  equation  8i 
is  satisfied  with  sufficient 

ACCUHACY. 

EPS2 

MAIN 

SMALL  PARAMETER  INTRODUCED  IN 
EQUATION  89,  WHICH  DETERMINES 
whether  Last  term  on  MiOHT- 
HaND  side  of  equation  13AIS 
significant. 

ERR 

MAIN.  MaSMON. 
OENSTy 

BIT(I)I  A TRUTh/FALSE  SWITCH 

TO  TELL  The  main  PROGRAM 

IF  Certain  convergence 

WAS  MET. 

S 

main,  prntal. 

MASMON 

gravitational  constant, 

OCCUhRING  in  equation  6. 

qmax 

MAIN.  MDhEN 

CUT-UFF  KaRAMTeR.  labeled  as 

"To"  IN  EQUATION  96,  FOR  SOLU- 
TION UF  EQUATION  13A. 
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mhere 

VARIABLES  USED 


DESCRIPTION 


IMAX 


MaIN«  PSAQS*  INITALt 
PRNTAL*  MASNONf  OENSTYf 
MOMEN 


NUMBER  nP  CELLS  INTO  MHICH 
CUMPUTAtIONAL  OnIO  is  OlVIOEO 
IN  X-DIRECTION, 


IMl 

main* 

INITAL 

IM2 

main* 

INITAL 

IM21 

main* 

INITAL 

IM22 

main* 

INITAL 

- IMAX  * I 
> IMAX  * 2 

■ 2*rMAX  ♦ I 

■ 2«IMAX  ♦ 2 


I PL 
ISL 
12 


main*  prntal*  counter 

MOMEN 

main*  prntal*  cuunter 

DENSTy 

Main*  inital*  masmon  ■ 2*imax 


UMAX 


main*  psaos*  inital* 
PRNTAL*  MaSMON*  OENSTY* 
MOMEN 


number  of  cells  into 

NMICH  COMPUTaTIUNAL 
(>RID  IS  DIVIUEO 
IN  Z-DIRECTION, 


JMl 

main*  INITAL 

m JMAX  * 1 

JM2 

main*  inital 

■ JMAX  * 2 

JN21 

main,  inital 

■ ?*JMAX  ♦ 1 

JM22 

main*  inital 

■ 2*JMAX  ♦ 2 

J2 

main*  inital*  masmon 

■ 2*JMAX 

M(20*40T 

main*  prntal* 

MASMON*  OENSTY 

MASS  IN  Each  CELL  Rij 

MOMX (20*401 

main,  prntal* 

MaSMUN*  MOMEN 

x-component  of  momentum 

IN  EACH  cell  Rij 

MOMZ(20*40) 

main*  prntal* 

MASMON*  MOMEN 

Z-CUMPONCNT  of  momentum 

IN  EACH  cell  Rij 

M2VX(20*40) 

main*  OENSTY* 

MOMEN 

CORRECTION  TO  X-COMRONENT  OF 
MOMENTUM  DUE  To  DENSITT  CONSTRAINT* 
labeled  "(Ap  1,"  IN  EQUATION  87A. 

M2VZ(20*40) 

Main*  oensty* 

MOMEN 

correction  to  2-comronent  of 
momentum  due  To  density  constraint* 

LABELED  " ( Ap  ),  ••  IN  EQUATION  87B. 

N 

main*  prntal 

COUNTER 

nmax 

main 

MAXIMUM  FOR  the  N COUNTER, 

ONE 

main*  MOMEN 

the  value  1.0 

P(20»40) 

Main,  psaos* 
oensty*  MOMEN 

coefficients  OfYlNO  EFFECT  OF 
DIFFUSION  IN  a.OIRECTIUN  OF  MASS 

AND  Z*C0MR0NENT  of  mUMENTUM* 

6IVEN  BY  EQUATIONS  eS  AND >6. 
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nhERE 

VARIABLES  USED 


MAiNt  PSAaS 


P1(20«A0)  PSAOS,  MOmEN 


O(AOtAO)  PSAQSt  OCNSTV 


OKAOtAOl  PSAQSt  HOMEN 


RHOiZOtAO) 


MAlNt  MASMONe 
OENSTy 

main*  OENSTY* 
MOMEN 

main*  PSAOS 


TM(A0*80)  MASMON 


TMOMX(A0*80i  MASMON*  MOMEN 


TMOMZ(A0*80)  MASMON*  MOMEN 


TPCON 


U(20*40) 


main*  PSAOS 


main*  PSAOS 


MAIN*  PSAOS* 
PBNTAL*  OENSTY 

main*  prntal* 

MaSMON*  MOMEN 


OtSCMIKTION 

V®»7»*  •hEBEAalS  STEP  SIZE  OF 
INOLPEnOENT  V'MIABUE  a 
IN  EQUATION  10A. 

COEFFICIENTS  GIVING  EFFECT  OF 
OIFFOSIUN  IN  X»DIRECTIUN  OF  X< 

component  of  momentum*  given  by 
EQUATIONS  98  AND  100. 

coefficients  Giving  effect  of 

OIFFOSIUN  IN  Z-OIPEdlON  OF  MASS 
AND  A^COHPONENT  OF  MOMENTUM* 
GIVEN  BY  EQUATIONS  69  AND  78. 

coefficients  giving  effect  of 

DIFFUSION  IN  Z.OIHECTIoN  OF  Z- 

component  of  momentum*  given  bT 
EQUATIONS  99  AND  101. 

DENSITY  IN  CELL  Rj]  • GIVEN  BY 
EQUATION  80. 

chapactemistic  density  of  fluid* 
labeled  "fs"  In  text. 


•>/A«*  PHEHEAa  IS  STEP  SIZE  OF 
independent  variable  a IN  EQUATION  10A. 

MASS  IN  Each  cell  of  eatenoeo 
ORID^  after  convection*  DENOTED 
BYm\j  In  EQUATION  53A. 

x-component  of  momentum  in  each 
cell  of  extended  grid  after  con- 
vection* denoted  BYh*x,jj*N  EQUATION53B. 

2-COMPONENT  OF  MOMENTUM  lY  EACH 
CELL  OF  EXTENDED  GRID  AFTER  CON- 
VECTION* denotfd  by  u*,  in  equation  B3C. 

2«/Aa/jr*  a mere  A a IS  STEP  SIZE  OF 
independent  variable  a IN  EQUATION  10A. 

2i/ST.  rmeRE  Aa  IS  STEP  SIZE  OF 
independent  variable  a IN  EQUATION  IDA. 

the  value  2.U 


X-COMPONENT  OF  VELOCITY  IN  CELL  Rji 
GIVEN  by  EQUATION  38A. 
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VARIABLES 

■HERE 

USED 

description 

V(?0««0) 

prntal.  oensty 

MOMEN 

* 

ouantity  mmich  describes  effect 
OF  density  constraint  UN 
momentum  given  by  equation  26. 

vmax 

main*  oensty* 

MOMEN 

maximum  UF  Vij  OVER  CELLS  Rn  * 
LABELED  ••  V+  ••  IN  EQUATION  91. 

VI (20**0) 

oensty*  MOMEN 

SCALED  values  OF  Vn  TU  INSURE 
stability  of  algorithm* 

GIVEN  BY  EQUATION  94. 

M(20**0> 

main*  prntal* 

MaSMON*  MOMEN 

Z-COMPONENT  of  VELOCITY  IN 

CELL  Rjj  GIVEN  by  EQUATION  38B. 

XMH(81) 

PSAOS*  inital* 
MASMON*  oensty 

X-COOHOIN*TES  OF  LEFT-HAND 

SIDES  OF  CELLS  IN  EXTENDED 

computational  grid*  given 

BY  EQUATIONS  35B  AND  43A. 

XPH(60) 

PSAOS*  inital* 

MASMON*  oensty 

X-C0UR01N«TES  oF  RIGHT-HAND 
SIDES  OF  CELLS  IN  EXTENDED 

computational  grid*  given 

BY  EQUATIONS  35B  AND  43A. 

ZMH(SI) 

PSAOS*  inital* 
MaSMON*  oensty 

prntal 

* z-coordtnates  of  bottoms 

OF  CELLS  IN  EXTENDED 

computational  grid* 

GIVEN  By  equations  35B 
AND  43B. 

zo 

main,  prntal* 
oensty*  MOMEN 

The  value  0. 

ZPH(BO) 

PSAOS*  INITAL* 
MASMON*  OENSTY 

PRNTAL 

* Z-COOROINATES  OF 

TOPS  OF  cells  in 

EATENOEO  COMPUTr- 

riONAL  6R10t  GIVEM 
BY  EQUATIONS  35B  AND  43B. 
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PAGE  1 

ROOCRSI  PROC  OPTIONSCMAIN) I /•  WATER  WAVES  6/25/78  •/  10 

20 

OCL  (ATANtSORTl  RUILTINl  30 

OCL  (INITAL*  PRNTAL*  PSAQSt  MASMONt  OENStY*  ROMENI  ENTHYI  AO 

OCL  OISKl  FILE  sequential  RCCOROI 

OCL  01SK2  file  sequential  MECOROI 


OCL  <Q«Q1)  (AOtAO)  float  BIN  EXTI 

OCL  (PtPl*  M*MOMX*MUNZtNHO«U«V»*t  VI  ) (20tA0)  FLUAT  BIN  EXTI 

OCL  (M2VXi  M2VZ) (20*401  FLOAT  BIN  EXTI 

OCL  (TMOMXt  TM0MZ«  TM>  140*80)  FLOAT  BIN  EXTI 

OCL  (DX(20)*  OZI40)*  XMH(81) *ZMH(8I) * XPh (60) *ZPH (60 ) I 

float  BtN  EXTI 

OCL  (OA*  00*  OT*  EPS*  EPSl*  EPS2*  8*  OHaX*  PCaN*  RHOC*  SODA* 
TPCONt  TSOA*  VMAX)  float  BIN  EXTI 
OCL  (zo*  ONE*  TWO*  PI)  Float  bin  exti 

OCL  (lMAX*JMAXElSL*lPL*NtIMl*IN2*ZM21tlM22*12,JMI*jN2*JM2I*JN22E 
J2)  FIXED  BINOl)  EXTl 
OCL  (I*J*NHAX)  fixed  BIN(31)| 

OCL  (DEBUG*  ERR)  BIT(|)  EXTI 


100 

no 

120 

130 

140 

ISO 

160 

170 

180 


ON  underflow I 190 

/••  200 

DEBUG  ■ 'PBI  210 

••/  220 

DEBUG  ■ •O'BI  230 

ERR  > *0*01  240 

ISL*  IPL  > 01  2S0 

ZO  ■ 0,01  260 

ONE  ■ I. 01  270 

TWO  > 2.01  280 

PI  * 4,0*ATAN(  ONE  ) I 290 


/•  —————  input  CASK  r/A/r« 

NMAX  ■ 201 

IMAX  ■ lOl 
JMAX  ■ 201 
NaOl 
GaUNCI 

EPS  a l.OE-SI 
fPSl  a l.nF-21 
EPS2  a 1.0t-.)l 


•/ 


340 

350 

360 

370 

380 
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OA  > O.OSEOI 

RHOC  ■1.01 

OMAX  ■ 10.01 

SODA  ■ SORT!  Da  l< 

TSUA  ■ 2.0  • SODA  I 

PCON  ■ SORT!  Oa/PI  >1 

TPCON  ■ TRO  • PCON I 


IMl  ■ IMAX  * II 
IM2  ■ IMaX  ♦ 21 
12  ■ 2*IMAXI 
IM22  ■ 12*21 
IM21  ■ 12  « II 
JMl  ■ JMaX  • It 
JM2  ■ JMaX  • 21 
J2  ■ 2*JMAXI 
JM21  ■ J2  * II 
JM22  ■ J2  • 21 


or  ■ O.OSEOI 


NONX*NOM2.M.RHO.U»R  ■ 201 


OX  ■ O.SEOl 
DZ  ■ O.SEOl 
00  1«1  TO  IMAXi 
00  Jal  TO  A I 
NdtJI  ■ 0.2SE0I 
ENOI  END I 

DO  1«1  TO  Al 
00  Jal  TO  lAI 
M(1*J)  a 0.2SE0I 
ENOI  ENOI 
00  lal  TO  IMAXI 
00  Jal  TO  4 1 

MOMXd.J)  a -lO.OaMdtJII 
ENOI  ENOI 
DO  lal  TO  41 
00  JaS  TO  141 
MOMZdtJ)  a alO.O*MdtJ)l 
ENOI  ENOI 
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CALL  INITALI 
CALL  PRNTALI 

/•  ...  ...  ...  ...  ...  ...  ...  .—  — — — ...  ...  .—  ...  ...  •/ 

CALL  PSAOSI 

/•••••••••••• 

IF  ONE  ■ 1.0  Then  so  to  FInII 

•••••••/ 

/•  ... ...  ...  ... •/ 

/•  ...  ...  ...  ...  ...  ...  ...  ...  ...  ...  ...  ...  ...  ...  ...  */ 

NEXTIHEl 

N a N*ll 

IF  N > NHAX  Then  go  to  FINII 
CALL  MaSHONI 

IF  ERR  then  60  TO  ERROuTi 
CALL  OENSTYI 

IF  ERR  Then  so  to  ERROuTI 
IF  WHAX  <■  EPS?  Then  001 
00  lal  TO  IHAXI 
00  jal  TO  JMAXI 

MOMX(ItJ)  a MOMX(I*Jl  « MZVXIItJ)  / OTI 
NOMZIItJ)  a MOMZdtJ)  • H2VZ(I,J)  / OTI 
ENOI  ENOI 
60  TO  OUT  I 

ENni 

CALL  MOMENI 


/a 

OUTt 

/a 


•/ 


730 

740 

750 

760 

770 

780 

790 

800 

610 

620 

630 

840 

850 

660 

670 

680 

890 

900 

910 

920 

930 

940 

950 

960 

970 

980 

990 

1000 

1010 

1020 

1030 

1040 

1050 

1060 

1070 


IF  MOO(NilO)aO 

THEN 

1080 

IF  MODINt  2la0 

Then 

1100 

IF  MOOtN*  5)ao 

THEN 

1090 

1110 

CALL  PRNTALI 

1120 
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PA0E  A 

IF  N»l  THEN  001 

WRITE  FILE(01SK2)  FROMIM)! 

WRITE  FILEIOISkZ)  FR0M(H0MX)I 
WRITE  FILE(D1SK2>  FHOmwOMZ)  » 

CLOSE  FILE (DISKS) I 
END  I 

IF  N<NMAX  THEN  GO  TO  NEXTIHEI 

WRITE  FILE<0ISK2)  FR0M(H)I 
WRITE  FILECOISkZI  FR0M(M0HX)I 
WRITE  FILE(OISk?)  FROM(HOHZ)t 
CLOSE  FILE(0ISK2)I 

60  TO  FINII 

CRROUTI 

PUT  SKIP  LISTdHELP'l  I 
CALL  PRNTALI 

FINll 

end  ROGERS! 

//GeSYSPRINT  DO  OUTLIMnSOOOO 

//G. DISKS  DD  0SN-RCP.FAV,H06D3e0ISPnULD* 

//  DCB«(RECFMRFtLRECL>3200«BLKSIZE>3S00) •SPAC£a(3S00t (6tS) eRLSE) 
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initali  phoci 

INIT 

10 

INIT 

20 

OCL  07(40)  t XMH(ai)  tZHH(81)  • 

XPH(60>  tZPH(<iO)  ) 

INIT 

30 

FLOAT  BIN  EXTI 

INIT 

♦ 0 

OCL  (lMAX«jMAXtTMl,lM2«lM21tIM22tIZ*JMltJM2.jH21«JM22t 

INIT 

so 

J2)  FIXEO  BINOl)  EXTI 

INIT 

60 

OCL  (It  J»  PIXEO  BIN(31 » 1 

INIT 

70 

INIT 

80 

XMH(l)sOi 

INIT 

90 

00  Izl  TO  IMAXI 

INIT 

100 

XMH(I«1)  s XMH(I)  « 0X(I)l 

INIT 

110 

ENDI 

INIT 

120 

00  I»IM2  TO  IM21I 

INIT 

130 

XMM(I)  ■ 24XMH(IMn  - XMH.IM22-I)» 

INIT 

UO 

endi 

INIT 

ISO 

00  I«l  TO  I2( 

INIT 

160 

XPH(I)  s XMH(I*1) t 

INIT 

170 

ENOI 

INIT 

180 

INIT 

190 

ZMH(1)>0I 

INIT 

200 

or  J»1  TO  jmaxi 

INIT 

210 

Zmh(j»i)  « ZMm(J)  * 0Z(J)l 

INIT 

220 

ENOI 

INIT 

230 

00  JbJM2  TO  JH21I 

INIT 

240 

ZMM(J)  « 2*ZNH(JM1)  - ZMH(JM22-J)I 

INIT 

2S0 

ENOI 

INIT 

260 

PUT  SKIP  LIST(i  j XMH  aHH 

Zmh  ZPH»)I 

INIT 

270 

00  Jsl  TO  J2I 

INIT 

280 

ZPM(J)  » ZMH( J*n  1 

INIT 

290 

PUT  SKIP  EDIT(J,  XMM(J)t  XPri(J)t  ZMM(J)f  2PH(J)) 

INIT 

300 

(F(4) t (4)F(intl) > 1 

INIT 

310 

ENOI 

INIT 

320 

INIT 

330 

ENO  INITaLI 

INIT 

340 

106 
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prntali  proci 

PRNT 

OCL  {M.MOHXtMOM7t  UtV,W  )(20f40)  FLOAT  BIN  CXil  PRNT 

DCL  (ZMHCBDt  ZPH<B0))  FLOAT  BIN  EXTl  PRNT 

OCL  (EPS#  G#  Two#  20  I FLOAT  BIN  tXT*  PRNT 

OCL  (IMAX#IPL*ISL«JMAX#N  ) FIXED  BlNOll  EXTl  PRNT 

OCL  (I*J#I1*I2)  FIXED  BIN#  (MIN)  BUlLTINI  PRNT 

OCL  (MT#  MIJ#  mOMXT#  M0M2T#  ET»  ETl#  ET2)  INIt(O)  FLOaT  BINI  PRNT 

DCL  SMJ(IMAX)  float  BIN!  PRNT 

SMJnZOI  RRNT 

PRNT 

11  • MIN(10#IMaX) I RRNT 

12  ■ MIN(IMAX#20) I RRNT 

RUT  page  EOIT(iFOH  N ■ • #N)  ( A#F (4) ) I PRNT 

RUT  EOIT(»#  •_LOOR  TOTAL  ■ •#ISL#'#  ♦_L00P  TOTAL  ■ •»IPL)  RRNT 

(A#F(B))«  RRNT 

00  I«1  TO  IMaXI  RRNT 

DO  J«l  To  JMAXI  RRnT 

MIJ  a M(I#J) » RRNT 

IF  MIJ<ERS  THEN  U(I#J)#  »I(I#J)  a Z0»  RRNT 

ELSE  001  RRNT 

U(I#J)  = MOMX(I«J)/MIj|  PRNT 

B(I#J)  ■ M0M2(I*J)/MIjl  PRNT 

END!  RRNT 

MT  a MT  ♦ MIJI  RRNT 

MOMXT  a MOMXT  « MOMX(I«J)f  RRNT 

MOMZT  a MOMZT  » HOMZ(I#J)l  RRNT 

RRNT 

ETl  a ETl  ♦ MIj  as*  (ZPH(J)*ZMH(J) )/T»OI  RRNT 

ET2  a ET2  ♦ Mlj  • (U(l.J)a*2  ♦ W ( I # Jl ••2) /TwO|  RRNT 

END I END  I RRNT 

RRNT 

ET  • ETl  ♦ ET2I  RRNT 

RRNT 

NEOZi  RRNT 

00  lal  TO  IMAXI  RRNT 

00  Jal  TO  JMAXI  RRNT 

SMJ(I)  a SMJd)  * M(1#J)I  RRNT 

ENUI  END I RRNT 

RRNT 

PUT  SKIP(2)  LISTC  M Ml  PRNT 

PUT  DATA(MT)I  RRNT 

DO  JaJHAX  TO  1 RY  >11  PRNT 

RUT  SKIP  EDIT(J,  (M(I#J)  DO  I«1  TO  Il))(H2).lO  F(12«5m  PRNT 

ENDI  PRNT 


10 

20 

30 

40 

SO 

60 

70 

80 

90 

100 

110 

120 

130 

140 

ISO 

160 

170 

180 

190 

200 

210 

220 

230 

240 

2S0 

260 

270 

260 

290 

300 

310 

320 

330 

340 

350 

360 

370 

380 

390 

400 

610 

420 

430 

440 

450 
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IF  IMAX>10  THEN  OOl 

PRNT 

A60 

PUT  SKIP! 

PRNT 

♦ 70 

DO  JaJMAX  TO  1 RY  >11 

PRNT 

♦80 

PUT  SKIP  EDIT! J« (H(ItJ)  00 

I«ll  TO 

12) ) (F|2>  ,10  F|l2.t>l  ) 1 

PRNT 

♦90 

ENOi 

PRNT 

SOO 

ENOI 

PRNT 

510 

PUT  SKIP<2)  EOIT(  SMJ  )(X(2ltlO  F(l2t5)l| 

PRNT 

520 

PUT  SK1PI2)  LIST(*  MOMX  t)! 

PRNT 

530 

PUT  OATA(MOMXT)» 

PRNT 

5^0 

00  JaJMAX  TO  1 RY  •!! 

PRNT 

550 

PUT  SKIP  E0lT(J.(H0MX(ItJ) 

UO  lal 

TO  ID  ) |F|2)  ,10  FlUtS)  ) 1 

PRNT 

560 

ENOI 

\ 

PRNT 

570 

IF  IMAX>10  THEN  001 

PRNT 

580 

PUT  SKIPI 

PRNT 

590 

DO  JaJMAX  TO  1 RY  -II 

PRNT 

600 

PUT  SKIP  EDIT(J,(M0MX  ItJ) 

UO  lall 

TO  12) ) IFI2) ,10  F 112,5) ) 1 

PRNT 

610 

ENOI 

PRNT 

620 

ENOI 

PRNT 

630 

PUT  SKIP(2)  LISTC  M0M2  'll 

PRNT 

660 

PUT  OATAINOMZTl  1 

PRNT 

650 

00  JaJMAX  TO  1 RY  -11 

PRNT 

660 

PUT  SKIP  EDITIJ,  (MOMZdtJ) 

00  lal 

TO  ID  ) IFI2)  ,10  FI12,5)  ) 1 

PRNT 

670 

ENOI 

PRNT 

680 

IF  IMAX>10  THEN  001 

PRNT 

690 

PUT  SKIPI 

PRNT 

700 

00  JaUMA*  TO  1 RY  -11 

PRNT 

710 

PUT  SKIP  EDITIJ, fMOMZlIfJ) 

00  lall 

TO  12) ) IFI2) ,10  Fll2,5) ) 1 

PRNT 

720 

ENOI 

PRNT 

730 

ENOI 

PRNT 

760 

PUT  SKIP(2)  OATAIFTlt  ETZt 

ET)  1 

PRNT 

750 

IF  Nao  Then  petujni 

PRNT 

760 

PUT  SKIP(2)  LIST(*  U») 

1 

PRNT 

770 

00  JaJMAX  TO  1 BY  -11 

PRNT 

780 

PUT  SKIP  EDITIJ,  (Ud.JI  DO 

lal  TO 

ID  ) IF  12)  ,lo  FI12,5)  ) 1 

PRNT 

790 

ENOI 

PRNT 

800 

IF  IMAX>10  then  001 

PRNT 

810 

PUT  SKIPI 

PRNT 

820 

00  JaJMAX  TO  1 RY  -11 

PRNT 

830 

PUT  SKIP  EDITIJ, (U(I»J)  00 

lall  TO 

12) ) IFI2) ,10  FI12,5) ) 1 

PRNT 

860 

ENOI 

PRNT 

850 

ENOI 

PRNT 

860 

PUT  SKIP(2»  LIST!'  rtd 1 

PRNT 

870 

DO  JaJMAX  TO  1 RY  -11 

PRNT 

880 

PUT  SKIP  EOITIJ,  |Wd,J)  DO 

lal  TO 

11)  ) IF  12) ,10  FI12,5) ) 1 

PRNT 

890 

ENOI 

PRNT 

900 

108 
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1 


i I 


! i 

i. 

a.  ‘ * 

I 

{ :* 


I- 


j. 


f ; 

I 


M 
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IF  IMAX>10  THEN  001 
PUT  SKIPI 

00  JnJMAX  to  1 BY  -II 

PUT  SKIP  EOIT(J,(wiItJ)  00  I»ll  TO  I2)I<F(2J»I0  F(12t5»IM 
ENOI 

FNOI 

PUT  SKIPI2)  LIST!'  V«) I 

DO  J«JMAX  TO  1 BY  -II 

PUT  SKIP  EDiTIJ.lVdtJ)  00  l»l  TO  UI>(F(2),lo  F(12ti))ll 
ENOI 

IF  IMAX>10  THEN  OOl 
PUT  SKIPI 

00  JnJMAX  TO  I BY  -II 

PUT  SKIP  EOITCJ,  (V(I.J)  DO  1 = 11  TO  12))(F(2)tiO  Fll2it>))l 


ENOI 

END  PPNTaLI 


ENOI 


PRNT  910 
PHNT  920 
PRNT  930 
PRNT  940 
PRNT  950 
PRNT  960 
PRNT  970 
PRNT  980 
PRNT  990 
PRNTIOOO 
PRNTlOlO 
PRNT1020 
PRNT1030 
PRNT1040 
PRNT1050 
PRNT1060 
PRNT1070 
PRNT1080 
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PSAQSI  PROCI 

OCL  (ERFCt  EXP  ) 
DCL  (OtOl) 

DCL  (PtPl 
DCL  (OX(20)f 


P*8E  9 

RUILTINJ 
(40, AO)  FLOAT  BIN  E*TI 

) (20,40)  FLOAT  BIN  EXT  i 
OZ(40lt  XMH(Sl) •ZMH(ai> , XPH{60) ,ZPH(bO) ) 


OCL  (IMAX.  JMAX)  FIXED  BInOI)  EXT) 
OCL  (PCONt  SQOa,  TPCON,  TSDA,  TXPm# 
OCL  PICON  FLOAT  PlNI 
DCL  (I,  Jt  K*  L)  FIXED  aiN|31)« 

PUT  PAGE  LISTI*  THE  P(I,J)"S*)« 
TXPM  » 2,0  • XPHI  IMAX  ) I 


DO  I-l  TO  IMAXI 
DO  K«1  TO  iMAXI 
IF  l^.K  then  1)01 

P(I,K)  « PFUN(  (XMH(K) ♦XMH ( I ) ) /TSOA 
♦ PFUNI  (XPH(M ♦XHH( 1 ) ) /TbUA 


FLOAT  BtN  EXT» 

TWO  ) float  bin  EXTI 


ENOI 


END! 

END! 


DO  I«1  TO  IMAXI 
P(ItI)>OX(I) I 

DO  K«1  TO  IMAXI 
IF  K-.»I  THEN  P(I*I) 
ENOI 
ENOI 


» P(If  I)  - P(ItK) I 


10 

20 

30 

40 

50 

60 

70 

80 

90 


psao 

PSAO 
PSaQ 
PSAO 
PSAO 
PSAO 
PSAO 
PSAO 
PSAO 
PSAO  100 
PSAO  110 
PSAO  120 
PSAO  130 
PSaQ  140 
PSAO  ISO 
PSAO  160 
PSAO  170 
PSaO  180 
PSAO  190 


♦ 

PFUN( 

ITXPH  - 

XPmIK)  - XPmII) )/TSUA  I 

PSAO 

200 

• 

PFUN( 

ITXPH  - 

XMhIK)  - XOflll)  )/TSUA  ) 

PSAO 

210 

* 

PFUNI 

ITXPH  - 

XPhIK)  - XMnll) 

)/TSOA  ) 

PSAO 

220 

♦ 

PFUNI 

ITXPH  - 

XMmIK)  - XMMlI) )/TSUA  ) 

PSAO 

230 

• 

PFUNI 

IXPHIK) 

*XMMin»/TSUA  ) 

- PFUNI 

|XMHIk)*XPhI1) )/TSDA 

) 

I 

PSaQ 

240 

PSAO 

250 

IF  K>1  THEN 

PII,K)  a P|I,K) 

PSAO 

260 

PFUNI 

IXMHIk) 

-XPMI I ) ) /TSUA  ) 

- PFUNI 

|XPHIK)-XPhI I) )/TS04 

) 

PSAO 

270 

m 

PFUNI 

IXMHIk) 

-XMMII) )/TSUA  ) 

* PFUNI 

(XPHIk)-XMnII) )/TS0A 

) 

1 

PSAO 

280 

PSAO 

290 

ELSE 

PII*K)  a PII,K) 

♦ 

PSaO 

300 

PFUNI 

IXHHII) 

-XPMIK)  )/TSL<A  ) 

- PFUNI 

IXPHID-XPHIR)  )/TS0A 

) 

PSAO 

310 

m 

PFUNI 

IXMHI  I) 

•XMHIK) )/TSUA  ) 

♦ PFUNI 

IXPHln-XMFilK)  )/TS04 

) 

1 

PSAO 

320 

PSaO  330 
PSAO  340 
PSAO  350 
PSAO  360 
PSAO  370 
PSAO  380 
PSAO  390 
PSAO  400 
PSAO  410 
PSAO  420 
PSAO  430 
PSAO  440 


1 

't  * I 
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00  I>1  TO  IRAKI 

PUT  SKIP(2)  EOITI  I ) (FIS)  ) I 

PUT  EOITI  (P(I.J)  DO  jNl  TO  IRAK)  )(F(U«5))I 

END  I 

PUT  SKIPIA)  LISTI*  THE  UII»J)"S«)I 
00  jRl  TO  URAXI 
00  LRl  TO  JMAXI 

IF  J»L  then  Q|J,U)«  UZIJ)  - 2*PCUN  ♦ 

2.0*PFUN(  DZ(J)/TS0A  ) 

♦ PFUNI  2RH|J)/S0DA  I - 2.0*PFUNI  ( ZMh (J) IJ I ) /TSDA  ) 

* PFUNI  ZPHI J)/SQ0A  ) I 


ELSE  DO  I 

OIJtL)  ■ PFUNI  (ZmH(L) «ZHH(J) ) /TSDA  ) - 

PFUNI  IZPMIL) ♦ZMHI J) I /TSOA  I-  PFON(  I ZRH <L I *ZPN (J) ) /TSOA  > 
♦ PFUNI  IZPhIL)*ZPHIJ>)/TSOA  )l 


IF  J<L  then  OIJ.L)  a Q(JiL)  * PFUNI  I ZMH (L) -ZPH I J) I /TSOA  ) 
- PFUNI  (ZPHILI-ZPmIJ) I/TSOr  I • PFUNI  I Zmh IL ) -ZMM ( J) » /TSOA  ) 
* PFUNI  |ZPh(L)>ZHHIJ) )/TSDA  )l 


ELSE  OIJ.L)  « ulJ.L)  .PFUNI IZMH(J)-ZPH(LI I/TSOA) 

- PFUNI  (ZPhiJI-ZPHILI  )/TS0A  )-  PFUM  (ZRH  I,|) -ZRM  (L ) ) /TSOA  ) 
♦ PFUNI  IZPH(J)»ZMHIL) )/TSUA  )l 
ENOI 

END  I 
ENOI 

DO  1*1  TO  JMAXI 

PUT  SKIPI2)  EOITI  I ) (FIS)  ) I 

PUT  EOITI  (QII.J)  DO  J»l  TO  JMAXI  l(F(12,5))l 
ENOI 


PUT  PAGE  LISTI*  THE  P1II*J)"S*II 

00  iBl  TO  IMAXI 

PICON  ■ OXII)  - TPCONI 


P1(1*I)  ■ TISO*PFUN(  OXiD/TbOA  ) ♦ PICON 

- PFUNI  XMH|I)/SOOA  ) - PFUNI  XPH(U/’SOOA  I 

♦ ThORPFUNI  IXMHII) .aPHII) )/TSDa  ) 

- PFUNI  IXPHIIMAXI-XPHID )/bUOA  ) 

♦ TRO*PFUN(  ITXPH  - aPHiD-XMHII)  )/TS»UA  ) 

- PFUNI  IXPhIIMAX)  - XMHII)  l/SClUA  )l 


PSaO  ASO 
PSAO  A60 
PSaO  A70 
PSaO  480 
PSaO  490 
PSaO  500 
PSaO  510 
PSAO  520 
PSAO  530 
PSaO  540 
PSAO  550 
PSAO  560 
PSAO  570 
PSaO  580 
PSaO  590 
PSaO  600 
PSaO  610 
PSAO  620 
PSaO  630 
PSAO  640 
PSAO  650 
PSaO  660 
PSAO  670 
PSAO  680 
PSaO  690 
PSAO  TOO 
PSAO  710 
PSAO  720 
PSAO  730 
PSaO  740 
PSaQ  750 
PSAO  760 
PSaQ  770 
PSaO  780 
PSaO  790 
PSaO  800 
PSAO  810 
PSaQ  820 
PSAO  830 
PSAO  840 
PSAO  850 
PSAO  860 
PSAO  870 
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00  K«1  TO  IMAXI 

IF  ih«k  then  noi 

IF  K>I  THEN  PICON  a PFUN  I (XMH (K) -XPH ( I ») /TSOa  ) 

- PFUN(  (XPH|K»-XPh<I) )/TS0A  ) - PFUNI  (XMH (K ) .XHH ( 1 ) ) /TSOA  ) 

♦ PFUN(  (XPH(K)-XMH(n )/TSOA  )l 

ELSE  PICON  a PFONI  IXMH ( 1 1 -XPH (K ) » /tSOA  ) 

. PFUNI  <XPHII).XPh(K) )/TSOA  ) - PFUNI  I XMH 1 1 ) nXMH (K ) I /TSOA  ) 

♦ PFUNI  IXPHID-XMHIK)  l/TSOA  )( 

PKItK)  ■ Two  • PICON  - PII.K)I 
ENOI 

ENOt 

ENOI 

00  lal  TO  IMAXi 

PUT  SKIPI2)  EOITI  I ) IFI5)  II 

PUT  EOITI  IPIII.JI  00  J«1  TO  IMAX)  >IF(12t5))| 

ENOI 

PUT  SKIPIAI  LISTI'  The  Ul  1 1 » JI "S» ) I 
00  Jal  TO  JMAXI 

QKJaJ)  a OZ(J)  - TPCON  * TaO*PFUNI  OZ(J)/TSOa  > 

- PFUNI  ZMHIJI/SOOA  I - PFUNI  ZPHIJI/SOOA  ) 

♦ Two  • PFUNI  IZMHIJ)  ♦ZPnIjn/TSUA  )l 
00  L*1  TO  JMAXI 

IF  J-aL  THEN  QIIJaL)  ■ U|JtL)  - TW0*PFUN | I 2MM |J ) rZhM |L ) ) /TSOA) 

♦ TW0«PFUN|  IZmhi J) rZPhILI I /TSOA  I 

♦ TWO*PFUNI  IZPH(J)*ZMhIL) I/TSUA  I 

- TWO*PFUNI  IZPHIJ)  *ZPmIL)  I/TS'iA  •)  I 

ENOI 

ENOI 

00  lal  TO  JMAXI 

PUT  SKIPI2)  EDIT  I I ) IF  15)  II 

PUT  EOITI  |01li,j)  00  J«1  TO  JMAX)  )IFI12f5))l 
ENOI 

OCL  OA  float  BIN  ExTl 

OCL  lOXM.OZM)  Un  FLOAT  HIM 

OCL  IXf  CPt  CZ.  CM)I1M&X)  float  bin* 

IZ<  EP*  EM,  EZ) IJMaa)  float  bin* 

10X2*  DZZ*  DM),  DM2*  01*  03*  CZM*  CZP*  LZm*  TOa)  FLOAT  BINI 


OCL 

11 

FI 

XEO  BINI31) 1 

00  I 

«1 

TO 

IMAXI 

Xlt) 

a 

IXPHII)«XMHII» )/Twui 

ENOI 

00  I 

>2 

TO 

IMAXI 

OXMI 

1) 

s 

XII)  - XII-1)  1 

ENOI 

PSaO  BBO 
PSaQ  B90 
PSaO  900 
PSAO  910 
PSaO  920 
PSAO  930 
PSAO  9A0 
PSaO  9S0 
PSAO  960 
PSaO  970 
PSaO  980 
PSaO  990 
PSAO 
PSAO 
PSAO 
PSAQ1030 
PSaQIOAO 
PSAO 
PSAO 
PSaOIOTO 
PSAOIOBO 
PSA01090 
PSaOUOO 
PSAOlllO 
PSaO 
PSAO 
PSAO 
PSAO 
PSAO 
PSAO 
PSAO 
PSA01190 
PSaOUOO 
PSaOUIO 
psaouzo 

PSAQU30 

PSaOUPO 

PSaOUSO 

PSAOU60 

PSaOUTO 

psaoubo 

PSAOUOO 

PSAOUOO 

PSAOUIO 

PSAOU20 

PSAOU30 

PSaQUAO 


1000 

1010 

1020 


lOSO 

1060 


120 

130 

lAO 

ISO 

160 

170 

lao 
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m 


at- 


00  vfal  TO  jm 
2(J)  • (ZPH( 
ENOI 

00  J«2  TO  JHAXI 
02M(J)  B 2(J 
ENOI 


PAOE  12 

TO  UMAX  I 

(2PH( J) *2MH( J) I/TbOI 


2<J)  • 2IJ*1)I 


OXM(l)  « ox (1)1  02M(1)  a 02(1 

0XM(1MaX*1)  a 0X(IMAX>I  D2H 

TDA  a TNO*OAI 
DO  lal  TO  IMAXI 
0X2  a 0X(I)«a2  / 12,01 
DMl  a DXMd)**?! 

0M2  a 0XMII*1)**?I 

IF  lal  Then  D1  a 0X(1)**2I 

Else  di  a oxc-d**?* 

IF  lalMAX  THEN  03  a 0X(IMAX)a*2l 
ELSE  D3  a 0X(I*l)**2» 

01  a 01/12.01  U3  a 03/12.01 


02(1)1 

D2M(jmax*1)  a OZ(jaAX)i 


Ul  a IJJ  a UJ/IC.U. 

CP(I)  a T0A/(DM2*D3-DX2*DXH(I)*0Xm(1*1)*(01-DX2)*0xm{1.1)/0XH(I) 
CM(l)  a TDA/(OM1.01-OX2»OX'^(I)*OXH(l*l)*(03-UX2)*DXM<i)/OXM(l*l) 
CZ(l)  a l.OEO  - CP(I)  - CH(l) « 

ENOI 


00  jal  TO  JHAXI 

022  a DZ(J)aa2/12.0l 
I a 02M(J)**2I 
i a OZM(J«l)«a?l 
IF  jal  THEN  01  a 02(1)**2I 

ELSE  01  a 02(J-1)»*2I 
IF  JajMAX  Then  03  a DZ ( JMAX ) aaz I 
else  03  a 0Z|jal)*a2l 

01  a 01/12.01  03  a 03/12.01 


EP(J) 

EH(J) 

E2(J) 

ENOI 


iC.OI  uj  a uj^ic.ui 

T0A/(0M2.03-UZ2*02M(J)*0AM(J*1) ♦ (D1-0Z2) •OAM(J*l)/ 

OZM(J) ) I 

TOA/ (nMl*01-0Z2«uZM ( J) ao/M ( J*l ) ♦ (03-0Z2) *nxM ( J) / 

DZN(J*1) I I 

l.OEO  - EP(J)  - tM(J) I 


pSAOiaso 

PSA01360 

PSAQliTO 

psAOiaao 

PSA01390 

psaqiaoo 

PSaQIAIO 

PSA01A20 

PSAQ1A30 

PSaOIAAO 

PSaOUSO 

PSA01460 

PSaOIATO 

PSAOlAeO 

PSA01A90 

psaoisoo 

PSAOlblO 

PSA01S20 

PSAQ1S30 

PSaOISaO 

psaoisso 

) IPSAQ1S60 

) IPSaQISTO 
PSaQI&BO 
PSA01S90 
PSA01600 
PSA01610 
PSAO1620 
PSA01630 
PSA016AO 
PSAOlbSO 
PSA01660 
PSAQ1670 
PSA01680 
PSAQlt>90 
PSAOWOO 
PSaOITIO 
PSAO1720 
PSAQ1730 
PSAO17A0 
PSAQ17S0 
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PAGE  »3 

CZM  ■ CZ<1)*CM(1) » 

IF  CZM  >■  0.  Then  001 

PUlt*»  ■ 0,1 
P(ltl)  ■ CZM  • 0X(1) « 

P(lt2>  « CP(1I  • OX(l)  I 
Pl(l»l)  ■ (CZ(l)-CM(l)  )•  t>X(l)l 
END! 

00  I«2  TO  IMAX-11 
If  CZ(1>  » 0.  then  001 
P(I.«) » PI (It*)  ■ 0.1 
P(ItI-!)t  PKTtl-I)  ■ CH(I)*OX(I)| 
Plltl)  t PI (Itl)  ■ C/<I)*OX(l) « 

p(iti*i)t  Pi(i.i*n  « cp(I)*oa(i)i 

END  I 

END  I 


EZM  > E2(l>*fM(i)  I 
IF  fZM  >■  0,  Then  OOI 

u(i«*>t  unit*)  ■ u.t 
O(ltl)  « EZM  * DZd)  ( 
uldtl)  « (EZd)>EHd))  * i)Zd)< 

Qdt2)t  unit?)  * FPd)*()7d)l 

END  I 


PSAQ1760 

PSaOITTO 

PSAOirao 

PSAQ1790 

PSAQ1800 

PSAOISIO 

PI  di2)  * p(lt2)  I PSAOIH20 

PSAQIG30 
PSaOIGAO 
PSAOloSO 
PSaOIBGO 
PSaQIBTO 
PSaOIBSO 
PSAOIB90 
PSA01900 
PSA019IO 
PSA01920 
PSA01930 
PSaOIVAO 
PSA019SO 
PSAQ1960 
PS*019T0 
PSAOI980 
PSA01990 
PSA02000 
PSA020I0 
PSAO2U20 
PSa02030 
PSAO20A0 
PSAO20S0 
PS&Q2U60 
PSa)J2070 
PSA020SO 
PSa02090 
PSA02100 
PSAQ2U0 
PSAO2120 
PSAa2130 
PSAO21A0 
PSA021SO 
PSA02160 
PSA02170 
PSA02180 
PSA02190 
PSA02200 


CZPaCZ(lMAX)  • CP(IMAX)( 

IF  CZP  >»  0,  Then  001 

Il»IMAX-n  PdMftXt*)*  PKIMAAt*)  a 0.1 

PdMAXtIDt  PKIMAXtll)  * CHdMAA)  * OXCIMaX)! 
P(lMAXtlMAX)  a CZP  * OXdMAX)! 

PKIHAXtlMAX)  a (CZ  ( IMAX) -CP  (IMAA)  ) *DA  (IMAX)  I 
FNOI 


PUT  Page  LIST(tNFa  THE  P(ItJ)"sdl 

00  lal  TO  IMAXI 

PUT  SK1P(2)  E01T(  I ) <E (5)  ) • 

PUT  EDIT!  (PCItJ)  00  Jal  TO  IhaX)  )(Fd2.5))l 
ENDI 

PUT  page  LIST(»NPa  The  PnnJ)''S')l 

DO  lal  TO  IMAXI 

PUT  SKIP(2)  EO|T(  I )(E(5)  )l 

PUT  E0IT(  (Pldt.))  DO  jai  TO  IMAX)  )(Fd2t5))t 
ENOi 
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00  JaZ  TO 

ir  CZIJ)  >•  0 THEN  OUl 

0(Jf*>  t Q1  IJ«*>  ■ 0.1 

0(JtJ-U.  OllJtJ-ll  ■ EM(J)*DZ<J)I 

OlJtJ).  OllJtJ)  • EZIJI  • 0Z(J)I 

■ EP(J>  • UZ(J)I 

ENOi 

Enoi 

IF  EZ(JH*x)  >■  0 Then  001 

0(JHAX,«),  UKJHaX.*)  ■ 01 

0<jH«X,JMAX-n t 01 (JHAXtJHAX>ll  > EM(jMAX)  • UZ(JHAX)I 
Q(JMAX,JMAXI t 01 ( JHAX. JMAX)  * EZ<JHAX|  • OZ(JHAX)l 
ENOI 

PUT  SKIP(4)  LISTCNEH  T'lE  OII.Jl''S*ll 

00  I»l  TO  JHAXI 

PUT  SKIP(Z)  EDIT(  I »(F(5)  >1 

PUT  EOIT(  (0(1, J)  00  J»1  TO  JMaX)  )(F(U, 51)1 

ENOI 

PUT  SKIP(4)  LIST(*NEH  the  01 (I»J)«S'I I 

00  1«1  TO  jmaxi 

PUT  SKIP(2)  ED1T(  I )(F(5)  )l 

PUT  E0IT(  (OKI.J)  DO  J*!  TO  JHAX)  )(F(i2,5))» 

ENOI 


PFUNj  PROC  ( Z ) RFTURNS(FL0AT  BIN) ( 

OCL  (Z,PZ)  float  BIN! 

PZ  » PCON  * EXP(-(Z**2))  - Z*SO0A  * ERFC(Z)I 
RETURN)  PZ  )l 
ENO  PFUNt 


END  PSAOSI 
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PSA022IO 

PSAQ2220 

PSA02230 

PSA02240 

PSAO2250 

PSAO2260 

PSA02270 

PSA02260 

PSAI32200 

PSA02300 

PSAO2310 

PSA02320 

PSA(}2330 

PSaOZJAO 

PSAO2350 

PSA02360 

PSA02370 

PSa02380 

PSAQ2390 

PSA02400 

PSA02410 

PSAU2420 

PSA02430 

PSAQ2440 

PSAQ24S0 

PSA02460 

PSA02470 

PSAO2480 

PSA02490 

PSAO2S00 

PSA02510 

PSA02S20 

PSA02b30 

PSA02540 

•/PSA02550 

PSA02S60 

PSA02b70 
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INI  PROCI 

MASM 

10 

masm 

20 

OCL  (H,MORXtMOMZ«AHOtU«  M > (20* 

AO) 

FLOAT  BIN  EXTI 

masm 

30 

OCL  (TMOMX,  TMOMZf  TM»  (AOtBO)  FLOAT 

BIN 

EXTI 

mash 

AO 

DCL  (OX(20)*  0Z(*0)*  XMfiCai*  *ZmH(81»  * 

XPH(60) *ZPH(eO)  ) 

masm 

SO 

FLOAT  BtN  EXTI 

masm 

60 

DCL  (ETAM(2:.  ETaP(2).  XIM(2)*  XIP<2)I 

FLOAT  8INI 

MASM 

70 

DCL  <DT.  EPS*  G 

) 

float  bin  EXTI 

masm 

80 

DCL  (IMAX*  JMAX*  12*  J2  ) fixed  BINOI)  EXT  I 

MASM 

90 

OCL  ERR  BlTd)  EXTl 

masm 

100 

masm 

110 

OCL  (FLOOR*  MAX*  MIN)  dOILTINl 

masm 

120 

DCL (Al *A2*60T*TFMP.TMAX  *UDT*XM*XMA A*XP*ZM*ZMAX*2P ) float  3 INI 

masm 

130 

DCL(I *IP*I1*I3.J*J3*K*km*l.LM*LP.NP.N1) 

FIXED  BIN(3I) I 

masm 

140 

masm 

ISO 

GOT  ■ G*DT| 

MASM 

160 

MOMZ  ■ MOMZ  - GOT  • Ml 

masm 

170 

masm 

ISO 

DO  I«1  TO  IMAXI 

masm 

190 

DO  Jal  TO  JMAXI 

masm 

200 

IF  MU*0)<EPS  THFN  U(I*  J)  *<«  (I  tv))  >01 

masm 

210 

ELSF  DOI 

masm 

220 

U(I*J)  a MOMXdtJ)  / Md*J>l 

masm 

230 

M(I*J)  a MUMZd*J)  / Md*j)| 

masm 

240 

ENOI 

mash 

2S0 

RMO(I*J)  a Md,J)  / (Oxd)*DZ(J)  ) 1 

masm 

260 

ENDI  EnOI 

MASM 

270 

masm 

280 

TH*  TMOMX*  TMOmZ  a 01 

MASM 

290 

J DO  lal  TO  IMAXI 

masm 

300 

DO  jal  TO  JMAXI 

masm 

310 

••  SUBROUTINE  1 •••••/ 

masm 

320 

ZM  a ZMH(J)  ♦ ad*J)»OT) 

HASH 

330 

ZP  a ZPH(J)  ♦ W(I*J)aOTI 

MASM 

340 

NPaOl  ETAM,  CTaP*  XIM*  XIP  a 01 

masm 

3S0 

ZMAX  a ZPH( JMAX) 1 

masm 

360 

masm 

370 

IF  ZP  >«  ZMAX  then  DOI 

masm 

380 

IF  ZM  < ZMAX  Then  OuI 

MASM 

390 

NPall  ETAMd>»ZMI  tTAPd) 

=zmax; 

masm 

400 

ENDI 

MASM 

410 

FNDI 

masm 

420 

MASM 

430 

ELSE  DOI  /•  NO^  FOR  ZP<ZmaX  */ 

MASM 

440 

IF  ZM  <r  -ZMAX  THEN  DOI 

MASM 

4S0 

IF  ZP  > -ZMAX  Then  DOI 

MASM 

460 

NPaJI  FTAMd)=ZMAXI  ETAPd) 

a 2*ZMAX  ♦ ZMI 

MASM 

470 

ENUt 

masm 

480 

ENDI 

MASM 

490 
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ir  HP  > xM  THEN  noi 

IPall  XIM(l)aXH| 
ENDI 

else  r>Oi 

IP  a 21  XlH(l) 

ENDI 

/#•#•#•#•••  SUBHOUTINC  3 

00  Hal  TO  IP! 

00  KMal  TO  l?l 

IF  XMH(KH)<aXIH(ll) 
END  I 

PUT  SKIP  L1ST(«ERR0N-1 


XIPlIlaXPI 

XlP(2)  a API  AlP(l) 

#•••••••/ 

Them  if  XlM|ll)<aXPM(KM) 
NO  KM»)|  ou  TO  ERROUTI 


mash  500 

MASM  510 
MASM  520 
MASM  530 
MASM  540 
MASM  550 
MASM  560 
MASM  570 
MASM  580 
MASM  590 
MASM  600 
MASM  610 
MASM  620 
MASM  630 
MASM  640 
MASM  650 
MASM  660 
MASM  670 
MASM  680 
MASM  690 
•/MASM  700 
MASM  710 
MASM  720 
MASM  730 
MASM  740 
MASM  750 
MASM  760 
MASM  770 
MA''v  SO 
MA  790 
MASM  800 
MASM  810 
MASM  820 
MASM  830 
MASM  840 
MASM  850 
MASM  860 
MASM  870 
MASM  880 
MASM  890 
MASM  900 
MASM  910 
MASM  920 
THEN  bO  TO  NEXTII  MASM  930 
MASM  940 
MASM  950 


PA6E  16 

ELSE  001  /•  N08  FOR  2M  > >ZMaX  •/ 

IF  ZP  > 0 Then  001 
IF  ZM»a0  THEN  001 

NPall  ETAMCIIaZMI  ETAPdiaZPI 
ENDI 

ELSE  001 

Npa2l  ETAP(l)aZPI 
ETAM(2)  a ZaZHAX  ♦ 2M| 

ETAP(2)  a 2«ZMaXI 

ENOI 

ENOI  /•  FOR  ZP>0  •/ 

ELSE  001 

NPall  ETAH(l)  a ZM  ♦ 2«ZMAX I 
ETAPin  a ZP  ♦ 2*Z«AXI 

ENDI 

ENOI 

ENOI  /*  FOR  ZP<ZmaX  •/ 


IF  NP>0  THEN  DO! 

/•••••a  SUBROUTINE  ? #••••#••#/ 

UOT  a U(ItJ)  • nTI 

XMAX  a XPHIIMAXIt  T^AX  a 2*AMAXI 

TEMP  a (XMH(I)  ♦ U0T)/THAXI 

XM  a TMAX  • (Temp  - FLOOflITtMPI ) I 

TEMP  a (XPH(l)«UnT)/TMAA| 

XP  a TMAX  * (Temp  - FLOOR(TtMPni 


a TmaAI 
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P*6E  17 

NEXTII  DO  KPnKM  to  I2«  M*SM  960 

IF  XMH(KP)<XIP(IU  Then  IF  X IP < 1 1 )<«XPH (KP)  THEN  bU  TO  NEXT2I  MASH  970 
END  I hash  980 

PUT  SKIP  LIST(iER«0K-2  yO  KP»>I  SO  TO  EHHOUTI  MASH  990 

NEXT2I  00  Nlal  TO  NPI  MASHIOOO 

00  LH»1  to  J^^  MASHIOIO 

IF  2HM«LH)<»ET4M(N1)  then  if  ETaH(NI)<«ZPh<lM»  then  60  TO  NEXT3IMASH1020 
ENOI  MASHI030 

PUT  SKIP(2)  DATAdtJtNPtlPtlZiNliETAMvETAPtZHHtZPH) t MASHIOAO 


PUT  SKIP  LIST  ( IFR90H-3  .yO  LH»)» 
NEXT3J  00  UP*LH  TO  J2» 

IF  ZMH(LP) <FT4P<nI)  Then  if 
ENOI 

PUT  SKIP  LIST(iFHR0H-4  no  L*"  M 


SO  TO  EPPOUTl 


Then  if  etap<nI)<«zpm<lpj 

NO  Lt"  M so  TO  ERNOUTI 


NEXT4I 


00  K«KM  TO  KPI 
DO  L»LM  TO  LRI 

A1  • MIN(XPH(K)  fXlPdl)  ) - iiAX(XHH(K)  tXIHdl)  ) I 
A2  « MIN(ZPH«L) .ETaP (N1 ) ) - max (ZMM(L) *ETAH(N1) ) I 

TM(K.L)  « TM(K,L)  ♦ HH(J(I*J)  * Al  • A2I 
TMOMX(KtL)  « TMO-tXd^.LI  * HHOd, J)*0<I,J)  * AI  * / 
TMOMZCKtLi  » TmOhZ(X,L)  * NHOd.JM»d»J>  • Al  » i 
ENOI  ENOI  /•  ENO  OF  K AND  L LOOPS  < 
ENOI  ENOI  /•  end  of  II  AND  N1  LOOPS  •/ 


ENOI  /• 
ENU  LOOP 1 1 


J LOOP 
END  OF 


/*  tNO  OF  nP>0  00 

• / 

I LOOP  •/ 


I3»I2*ll  J3«J?*1I 


HASM1U50 

MASM1U60 

THEN  GO  TO  NEXT4IMASH1070 
HASH! 080 
MASH1090 
MASHIIOO 
MASHIUO 
MASH1120 
MASM1I30 
MASHIUO 
MASHIISO 

I MASHIUO 

MASHI170 
MASM1I80 

A^l  MASHII90 

A2I  MASH1200 

•/  MASMI2I0 

MASHI220 
MASMI230 
MASHI240 
MASMI250 
MASMI260 
MASHI270 
MASHI280 
MASM1290 
MASHI300 
MASHI3I0 

t3-I*J3-J)l  MASHI320 


00  I«I  TO  IHaXi  MASHI300 

00  0>I  TO  JMAXI  MASHI3I0 

MdfJ)  ■ TMd.J)  ♦ TMdJ-I,J)  * THd.JJ-J)  ♦ TMd3-I*  J3-J)  I MASHI320 

MOHXd, J1»TM0MX  d t J> -TmuNX (13-1 fO) ♦ThOHX (It J3-J) -TNUHa (13-1 , J3-J) | MASH 1 330 
MOHZ  d»J) »ThOH7 ( 1 1 J) ♦TMuMZ  d3-I» J)-TM0MZ ( 1 « J3-J)-TH0HX( I3-I t J3-J) IMASHI340 
ENOI  ENOI  MASHI3S0 

MASHI360 

PETUMNI  MASHI370 

EPPOUTI  EPHn'I'BI  MASHI380 

END  HASHONI  MASH1390 
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OENSTVI 

PROCI 

DENS 

10 

OCL 

(DEBUG*  ERR)  BlT(l)  EXT! 

OENS 

20 

OENS 

30 

OCL 

0(40*40)  float  bin  EXTI 

DENS 

40 

OCL 

(M*  M2VX*  M2VZ*  P*  RHO*  V*  Vl)  (20*40)  FLOAT  BIN 

EXTI  OENS 

so 

OCL 

(DX(20)*  DZ(40)*  XMH(ai) tZMH(81> • XPH (60 ) ,ZPH (GO ) ) 

OENS 

60 

float  bin  EXTI 

OENS 

70 

OCL 

(OA*  06*  EPSl*  RHOC*  VHAX*  ZO*  ONE*  TnO)  FLOaT  GIN  EXTI 

OENS 

80 

OENS 

90 

OCL 

(MAX*  MINI  RUILTINI 

OENS 

100 

OCL 

(IMAX*  UMAX*  ISL)  FIXED  BINI31)  EXTI 

DENS 

110 

OCL 

(I*  ICSTl*  J,  K*  L)  FIXED  9INI31)I 

OENS 

120 

OCL 

(PRHO*  PSUMl. PSUM2*  QRHU«0SUM1 *LSUH2*RH0MaX*XMPI *4MPJ 

OENS 

130 

) float  bini 

OENS 

140 

OENS 

150 

V*  M2VX*  M2VZ  ■ 01 

OENS 

160 

ICSTl  » 01 

OENS 

170 

VMAX 

« ZOI 

OENS 

180 

/• 

— •/OENS 

190 

STAR.LOOPi 

' ,,  IF  DEBUG  THEN  PUT  SKIP  LIST ( *ENTEPEU  STAR_LOOPi)| 

: f • ISL  • ISL  ♦ II 

II.  ICSTl  « ICSTl  ♦ lit 

DO  I«1  TO  IMAKI 
5 , 00  TO  UMAX  I 

' 1 RHOtl.J)  » M(I,J)  / tDX(I)*02<j»)l 

{ I IF  l«j«2  Then  rhohax  > HHOdf i>-rholi 

i * ELSE  RHOMAX  > MftXtHHOHAXf (KHOII.J)-RmOCI ) I 

ENOI  ENOI 


DENS  200 
OENS  210 
OENS  220 
OENS  230 
OENS  2A0 
DENS  2S0 
OENS  260 
OENS  270 
OENS  280 
OENS  290 


OENS  300 
OENS  310 
OENS  320 
OENS  330 
OENS  340 
OENS  350 
OENS  360 
DENS  370 
OENS  380 
DENS  390 
DENS  400 
OENS  410 
OENS  420 
OENS  430 
OENS  440 
OENS  4S0 
DENS  460 
DENS  470 
OENS  480 


IF  DEBUG  THEN  PUT  SKIP  DATA (ICSTl*  MmOmax)! 
IF  RHOMAX  <«EPS1  Then  GU  to  FINAL_MUMXA2I 


VHAX  > ZOI 
00  I>1  TO  IMAXI 
XMPI  • XMH(I)  * XPH(I) I 
DO  J-1  TO  UMAXI 
PSUHl*  PSUM2  « ZOI 
00  K«1  TO  iMAXi 

PRHO  ■ PIK*T>  • MAX(HH0|K,J)-RH0C»  ZU)  I 
PSUMl  ■ PSUMl  ♦ PRHOl 

PSUM2  « PSUM2  ♦ PrtHO*(XMPI  - XMh(K)  - XPH(K)I  / iHUl 
ENOI 

V(1*J)  > V(I«J)  « OA  • MAX(RHO(I«J>>HHUCt  ZO)  I 
VMAX  s MAX(V(I,J>t  VMAXM 

M(I.J)  « OZ(J)*(r)X(I)*MlN(RHO(I*J)  >HH0C)  ♦ PSUMDI 
N2VX(I*J)  > M2yXll*J)  * 0Z(U)*PSUM2I 
ENOI  ENOI 
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00  I«1  TO  IMaXI 
00  jNl  TO  JHAXI 

RHOdtJ)  * M(ItJ)  / (0X(IJ»02<J)  ) I 
EMOI  ENOI 

00  Jal  TO  JHAXI 

2HPJ  * ZMH(J)  ♦ /PH( J) I 

00  I>1  TO  IHAXt 

USUMlt  OSUM?  > 201 
00  L*1  TO  JMAXI 

OWHO  « OlL*J)*MAA{RMU(IiL)-KMOCt  20)1 
QSUHl  B QSUMl  * URHOI 

0SUM2  B USUH?  ♦ (jRHU*(2HPj-2MH(L>-2Pn«Ln/T*(OI 
END! 

M<IfJ»  B OXlI>*(D7(J>*MlN(RnOII»J) . KhOC)  ♦ USUM] ) I 
H2V2II«J)  a H2V2(I»J)  • 0X(1)*USUM2I 
ENOI  ENOI 

IF  ICSTl  > 999  THEN  D0« 

PUT  page  LIST(*rcSTl  > 999  "STOP"  •)! 

ERNBipBl  RETUHNI 
ENOI 

GO  TO  STaR.LOOPI 
FINAI MOMXAZi 

IF  DEBUG  Then  put  skip  list  I »EnTEHEO  FINaL_M0MXa2»I  I 
IF  DEBUG  then  PUT  SKIP  UATA(VMAX)I 
END  OENSTYI 


OENS  A90 
OENS  SOO 
OENS  510 
OENS  520 
DENS  530 
DENS  5A0 
DENS  550 
OENS  560 
OENS  5T0 
OENS  5B0 
DENS  590 
OENS  600 
OENS  610 
OENS  620 
OENS  630 
DENS  6A0 
OENS  650 
OENS  660 
OENS  670 
OENS  680 
OENS  690 
DENS  700 
OENS  710 
OENS  720 
OENS  730 
OENS  740 
OENS  750 
OENS  760 
OENS  770 
OENS  780 


120 


THE  JOHNS  HQMdNS  UNIVERSITY 

APPUED  PHYSICS  LABORATORY 

LAURCt  MARVLAMO 


PAGE  20 


MOMCNI  PROCI 

NOME 

10 

HOME 

20 

OCL  (EXP  laUILTINI 

HOME 

30 

OCL  (DEBUG*  ERR)  RIT(11  EXTI 

MOME 

40 

OCL  (0*01)  (40.40)  float  BIN  EXTI 

HOME 

so 

OCL  (PtPi*  MOMX«MOMZ«  UfVttif  VI  ) (20*40) 

float  bin 

EXTI  MOME 

60 

OCL  (H2VX*  N2VZ) (20*40)  FLOAT  BIN  EATI 

MOME 

70 

OCL  (TMOMX*  TMOMZ  ) (40*110)  FLOaT  BIN  EXTI 

MOME 

80 

OCL  (OX(20)*  DZ(40)f  DA*  06*  DT*  GMAX*  HHOC*  VMAX* 

ZO*  ONE) 

HOME 

90 

Float  btn  exti 

MOME 

100 

OCL  (IMAX*  JHAX*  IPL  ) FIXED  BlN(il)  EXTl 

MOME 

no 

MOME 

120 

OCL  (A*EP0G*ETFMP*PSUM1*QSUMI*VMAXU  ) float  BINI 

MOME 

130 

OCL  (I*  IP*  J>  K *L>  fixed  blNOni 

MOME 

140 

HOME 

150 

VMAXll  a 1.1  • VMAXI 

MOME 

160 

IPaOl 

MOME 

170 

A a 1.1  • VHAX/RHOCI 

MOME 

180 

DG  a OA/Al 

MOME 

190 

IF  DEBUG  THEN  PUT  SKIP  DATA(VMaX11«  A*  OG.  GMaX) I 

MOME 

200 

MOME 

210 

DO  lal  TO  IMAXI 

MOME 

220 

DO  Jal  TO  JMAXI 

MOME 

230 

VKl.J)  a V(I*J)  / VHAXllI 

MOME 

240 

TMUMXd.J)  a MOMXd.J)  * M2VXd*J)  / OTi 

MOME 

250 

TMOMZd.J)  a MOMZd.J)  ♦ HpvZd.J)  / DTI 

MOME 

260 

HOMXd.J).  MOMZd.J)  a ZOI 

MOME 

270 

ENOI  ENOI 

MOME 

280 

IF  DEBUG  Then  noi 

MOME 

290 

PUT  SKIP(2)  LISTC  TM0mX»)I 

MOME 

300 

00  JalO  TO  1 BY  -1  1 

MOME 

310 

PUT  SKIP  EDIT(J,  (TMUMXd.J)  DO  lal  TO  1U))(F(2)*10 

E (12*4) ) 1 

MOME 

320 

ENOI 

MOME 

330 

PUT  SKIP(2)  LlSTd  TM0MZ*)l 

MOME 

340 

DO  JalO  TO  1 BY  -11 

MOME 

350 

PUT  SKIP  E0IT(J*  (TMOMZd.J)  DO  lal  TO  10))(F(2>,10 

E (12*4) ) 1 

MOME 

360 

ENOI 

MOME 

370 

ENOI 

MOME 

380 

MOME 

390 

PLUS.LOOPi 

HOME 

400 

”lF  DEBUG  Then  put  SKIP  LIST ( 'ENTEREU  HLuS_LOOP')l 

MOME 

410 

IF  DEBUG  then  PUT  SKIP  i)ATA(IP)I 

MOME 

420 

IPL  a IPL  ♦ 11 

HOME 

430 

IF  (IP*1)»DG  > GMAX  THEN  DUl 

MOME 

440 

EPDG  a EXr(-iP*0G)l 

MOME 

450 

DO  lal  TO  iMAXI 

MOME 

460 

DO  Jal  TO  JMAXI 

MOME 

470 

MOMX(I.J)  a MOMXd.J)  ♦ TMOMX  (1  . J)  *EPOG  1 

MOME 

480 

MOMZd.J)  a MOMZd.J)  ♦ TMOMZ(I.J)*tPOGI 

MOME 

490 

ENOI  ENOI 

MOME 

bOO 

HETURNI 

MOME 

510 
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ENOt 


ETEMP  3 EXP(-IP*nG)*(ONE  - EXP(>06))I 
00  I3l  TO  iMAXt 
00  J3l  TO  JMAXI 

M0MZ(I«J)  » MrtMZdtJ)  ♦ THOMZ(I*«J)*ETfmP| 

MOMX(ItJ)  3 MOMX(l*JI  ♦ TMOHX ( 1 1 Jl RETEMPI 
ENOt  ENOI 

Ut  M 3 ZOI 
DO  l3l  TO  IMAXI 
DO  J3l  TO  JMAXI 
00  L»1  TO  JMAXI 

U(I*J)  3 U(I,J)  * Q(LtJ)*Vl (I*L)*TMOMX(ItL)/(OX(I)*OZ(L) ) » 
W(I.J)  3 KltJ)  ♦ Ql  (Lt  J»*Vl  (I*L»*TM0MZ(  It|  )/<DA(I)*OZ(L)  ) I 
ENOI 

ENOI  ENOI 

DO  l3l  TO  IMAXI 
DO  J3l  TO  JMAXI 
OSUMl,  PSUMl  3 ZOI 
00  K3l  TO  IMAXI 

PSUMl  3 PSUMl  * PI (K, I ) *0(Kf J) I 
QSUMl  3 QSUMl  4 P(XtI)*R(K<J)  I 
ENOI 

TMOMXdtJ)  3 (ONE-Vld«J))  * TMOMX(lfJ)  4 PSUMII- 

TMOMZdtJ)  = (0NE-Vld4J))  * TMOMZdfJI  4 OSUMII 
ENOI  ENOI 

IP  3 IP  4 II 
GO  TO  PLUS_LOOPI 
END  mOMENI 


MOME  S20 
MOME  S30 
MOME  !>A0 
MOME  550 
MOME  560 
MOME  570 
MOME  580 
MOME  590 
MOME  600 
MOME  610 
MOME  620 
MOME  630 
MOME  6A0 
MOME  650 
MOME  660 
MOME  670 
MOME  680 
MOME  690 
MOME  700 
MOME  710 
MOME  720 
MOME  730 
MOME  7A0 
MOME  750 
MOME  760 
MOME  770 
MOME  780 
MOME  790 
MOME  600 
MOME  810 
MOME  820 
MOME  830 
MOME  840 
MOME  8S0 
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0 


D 

D 

c 

gO 

+ 

f(*) 

g 

1 

J 

I 

J 

k 

I 


GLOSSARY 


Coefficient  In  representation  of  S^(*)  finite  difference 
operator 

Coefficient  In  representation  of  S^(*)  by  finite  difference 
operator 

Coefficient  In  representation  of  S^(*)  by  finite  difference 
operator 

Flow  domain,  domain  In  Eq.  33 
Computational  domain 


Coefficient  In  representation  of  S (•)  by  finite  difference 
operator 

Coefficient  In  representation  of  S (•)  by  finite  difference 
operator 

Coefficient  In  representation  of  S (•)  by  finite  difference 
operator 

Function  defined  In  Eq.  10b. 

Gravitational  constant 
Non*-negatlve  Integer 
Non-negative  Integer 

Number  of  computational  cells  In  horizontal  direction 
Number  of  computational  cells  In  vertical  direction 
Non-negative  Integer 
Non-negative  Integer 
Mass  In 

Intermediate  mass  In  R. . 
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“ij 

n 

P(*) 

Ik 


*11 

S(-) 

?(•) 


S (•) 

X 


S^(-) 


t 

u 


Mass  in  cell  at  end  of  time  step 

Non-negative  integer,  unit  vector  normal  to  boundary 
Function  defined  by  Eq.  73b 

Approximation  to  operator  for  horizontal  diffusion  with 
reflection  at  x ■ 0 and  x ■ X 

Approximation  to  operator  for  horizontal  diffusion  In 
Infinite  space 

Approximation  to  operator  for  vertical  diffusion  with  re- 
flection at  z 0 

Approximation  to  operator  for  vertical  diffusion  in  Infinite 
space 

Rectangle  defined  by  Eq.  36 

Operator  that  transforms  solution  at  one  time  Into  solution 
at  later  time 

Approximation  to  operator  that  updates  solution  by  one  time 
step 

Operator  for  diffusion  In  horizontal  direction  in  Infinite 
space 

Operator  for  diffusion  In  vertical  direction  in  infinite 
space 

Time 

Velocity  components,  horizontal  velocity  conqronent 


u 


u 


IJ 


Initial  velocity 

Intermediate  velocity 
Velocity  at  end  of  time  step 
Horizontal  component  of  velocity  in  R. 


IJ 


Quantity  defined  In  Eq.  12 
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Approslaatc  value  of  v In  cell 

Vertical  velocity  component 
Vertical  component  of  velocity  In 

Spatial  coordinates,  horizontal  spatial  coordinate 
Width  of  computational  domain 

Horizontal  coordinate  of  centers  of  computational  cells  In 
1th  column 

Horizontal  coordinate  of  right-hand  sides  of  computational 
cells  In  1th  column 

Widths  of  computational  cells  In  1th  column 

Vertical  spatial  coordinate 
Height  of  computational  domain 

Vertical  coordinate  of  centers  of  computational  cells  In 
jth  row 

Vertical  coordinate  of  tops  of  conq>utatlonal  cells  In  Jth 
row 

Heights  of  computational  cells  In  jth  row 

Pseudo-time  variable 

Step  of  pseudo-time  variable  a 

Pseudo-time  variable 

Cut-off  parameter  for  solution  of  Eq.  13a 

Stap  of  pseudo-tlsM  variable  y 
Laplaclan  operator 

Cut-off  to  keep  from  dividing  by  zero  In  Eq.  38 

Cut-off  In  Eq.  81  for  approximate  satisfaction  of  density 
constraint 


THE  JOHNS  HQMINS  UNIVENSITV 

AP(>UEO  PHYSICS  LABORATORY 

UUIMl  MORVUMO 


13«  is  solved 


Horizontal  component  of  momentum  In  R. 


Intermediate  horizontal  component  of  momentum  in  R. 


Vertical  component  of  momentum  in  R. 


Intermediate  vertical  component  of  momentum  In  R. 


Momentum  In  cell  R 


Density 


Density  of  liquid  phase 


Initial  density 


Intermediate  density 


Density  at  end  of  time  step 


Gradient 
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1 
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1 
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1 
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1 
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1 
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1 
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1 
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I 
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Lib. 

1 
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1 
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1 
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1 
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1 
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